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Chapter  1 

Introduction 


Numerical  models  of  physical  phenomena  often  employ  a  representation  of  the  computational  space  which  is  discrete. 
For  example,  a  differential  equation  may  be  solved  on  a  rectangular  grid  of  discrete  points  through  finite  difference 
methods.  While  these  representations  of  solutions  are  valuable,  it  is  often  desirable  to  approximate  the  solution  in  the 
regions  between  solution  points  or  "nodes”. 

For  slowly  varying  functions,  a  linear  interpolation  of  values  is  often  employed.  This  approach  may  provide  sufficient 
accuracy  in  some  cases.  However,  simulations  which  require  many  iterations  are  often  adversely  affected  by  the  error 
introduced  by  linear  interpolation.  For  example,  small  surface  faceting  errors  for  optical  surfaces,  introduced  early  in 
propagation  are  often  manifested  as  significant  aberrations. 

Our  current  application  for  a  generalized  spline  implementation  is  rooted  in  the  need  to  solve  coupled  physical  prob¬ 
lems  through  differing  methods.  For  example,  the  motivation  for  this  technical  report  was  the  need  to  combine  a 
heat  transfer  model  with  two  different  laser  beam  propagation  codes.  In  both  cases  the  temperature  distribution  in  the 
material  induces  a  change  in  the  optical  properties,  such  as  refractive  index.  The  coupled  code  for  the  laser  beam  prop¬ 
agation  simulation  require  values  of  the  optical  properties  on  either  a  differing  grid  spacing  (for  fast  Hankel  transform 
methods),  or  at  arbitrary  points  within  the  material  (for  ray  tracing  methods). 

We  have  proposed  a  modified  Monte -Carlo  approach  to  the  solution  of  the  radiative  transport  equation  which  has 
the  unique  feature  of  incorporating  refractive  index  gradients  within  a  multi-layer  biological  tissue  model.  In  the 
approach,  photon  trajectories  are  computed  using  a  solution  of  the  Eikonal  equation  (ray-tracing  methods)  rather  than 
linear  trajectories.  The  method  can  be  applied  to  the  specific  problem  of  incorporating  thermal  lensing  and  other 
non-linear  effects  in  turbid  media  (biological  tissues)  by  coupling  the  radiative  transport  solution  into  heat-transfer  and 
damage  models.  In  turn,  the  method  can  be  applied  in  the  establishment  of  laser  exposure  limits  for  tissue-penetrating 
wavelengths,  as  well  as  a  number  of  additional  applications  in  imaging  and  spectroscopy  as  well  as  vision  science. 
Presented  here  is  a  short  summary  of  the  theory  and  methods  for  the  implementation  of  a  spline  interpolation  suitable 
for  accurate  one  and  two-dimensional  functional  distributions.  Included  is  source  code  for  both  the  MatLab  and  C++ 
programming  languages.  Example  data  are  presented,  along  with  a  short  stability  and  error  analysis  for  problems  of 
recent  interest  to  our  research. 
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Chapter  2 


Definitions 


2.1  Splines 

Suppose  [a,b]  is  a  closed,  finite  interval.  We  select  N  distinct  points,  or  knots,  x  i , . . . ,  xN  strictly  between  a  and  b. 
That  is,  we  have 

a  =  x o  <  Xi  <  X2  <  ...  <  Xn  <  X,y+ 1  =  b.  (2.1) 

The  intervals  between  each  consecutive  pair,  x,  and  x,+i  for  /'  =  ()...  N  form  a  partition  of  [a,  b],  as  long  as  we  carefully 
define 


7/  =  [x;,X1+i) 


for  all  i  -0...N  -  1,  and 


In  -  [xjv,  Xiv+i] 

so  that  the  right  endpoint  b  is  in  the  partition,  which  we  call  A. 

Let  Pm  be  the  space  of  polynomials  of  order  m.  Recall  that 

Vm  =  jp(x)  =  ^  c,x'~ 1  d  e  r|  .  (2.2) 

P-i ,  for  example,  is  the  set  of  quadratic  functions. 

Finally,  let  m  be  an  integer,  -1  <  m  <  m  -  2.  This  will  determine  the  smoothness  we  would  like  at  our  knots.  If 
m  =  —  1,  then  we  do  not  require  the  spline,  or  any  of  its  derivatives,  to  be  continuous  at  the  knots. 

SCPm ;  m;  A)  denotes  the  space  of  polynomial  splines  of  order  m  with  knots  xn-  We  say  that  s  is  in  this  space  if 

there  exist  polynomials  So, . . . ,  Sn  e  Pm  such  that 

1.  s(x)  -  Sj(x )  Vx  e  /,,  /  =  A',  and 

2.  D’ Sj_\ (Xj)  =  D' Si(Xj)  for  j  =  0, 1, ... , m  and  for  i  —  1, . . . , N. 

Theorem  2.1.1.  [4]  S(Pm  ;  m;  A)  is  a  linear  (vector)  space  of  dimension  m  +  (m  -  m  -  1  )N. 

Definition  2.1.2.  [4]Let  the  partition  of[a,b]  and  m  be  given.  A  =  is  the  extended  partition 

associated  with  S(P,„;  m;  A)  provided  the  following  conditions  hold: 

1-  y  1  <  yi  <  •  •  •  <  y2m+(m—m—l)N 

2.  yi  <■■■  <  ym  <  a 

3.  b  <  y m+{m-m-\)N  +1  -  ’  "  -  y2m+(m-m- 1  )N 
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Ill  —  111  —  1 


III  —  III  —  1 


4-  ym+ 1  —  ‘  ‘  ‘  —  y ~ 


X\9...,Xi 


XN,  . . . , XN 


2.2  B-splines 


For  numerical  applications,  a  local,  symmetric  basis  is  more  suitable  than,  for  example,  a  one-sided  basis  [4],  We 
therefore  choose  to  use  a  B-spline.  Given  integers  i  and  m  >  0,  for  all  real  x  the  mth-order  B-spline  associated  with 
knots  yi  through  yI+m  is  given  by 


(- 1 T  tfi ,  •  ■  • ,  yi+m ]  (X  -  1  if  y;  <  y,+m  and  x  >  y 


0  otherwise 

where  [  ]  denotes  the  divided  difference.  For  a  function  /,  the  divided  difference  [t \ , . . .  ,tr+  i  ]/  is  given  by 


(2.3) 


det 


[tu...,tr+i]f 


1  U 

1  h 
1  : 


1  fr+l 


det 


l  1 
1 

1 


fi 

h 


1  tr+ 1 


f'r  /(g) 

/r1  /(f2) 


V+l 


//+ 1) 


‘1 

,/-l 


(2.4) 


fr-\  fr 
V+l  V+l 


Depending  on  spacing  of  the  knots,  however,  B-splines  can  be  very  large  or  very  small  which  can  yield  unfavorable 
results  computationally.  The  normalized  B-spline  associated  with  knots  y,, . . .  ,y;+m  is  given  by 


N'i'(x)  =  (yi+m  -y,-)GT(x). 


(2.5) 
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Chapter  3 


Approximation  in  one  dimension 


3.1  The  objective  function  in  one  dimension 

Once  we  have  determined  the  spline  basis  we  will  use,  we  need  to  determine  the  appropriate  coefficients  to  approximate 
the  function  n  at  arbitrary  points  in  the  (r,z)  plane. 

For  the  time  being,  we  will  consider  our  basis  for  r.  Let  the  i,h  basis  function  be  given  by  (f.  and  let  the  number  of 
basis  elements  (the  dimension)  be  K  =  m  +  (m  -  m  -  1  )N.  Then  our  estimate  h  for  the  function  n  at  an  arbitrary  point 
r  is  given  by 


K 

n(r)  =  ^  c,4>,(r)  (3.1) 

i=i 

and  c  =  (ci,C2, . . . ,  ck)  are  coefficients  to  be  determined. 

We  would  like  to  choose  these  coefficients  so  that  they  give  us  the  best  possible  approximation  to  the  true  values  of  the 
function.  Our  idea  of  what  kind  of  approximation  is  “best”  depends  on  the  particular  application.  However,  certainly 
we  would  like  the  approximation  to  come  close  to  the  known  values  n  at  the  gridpoints.  Secondly,  we  would  like  the 
approximation  to  remain  close  to  the  data  we  have  on  the  grid,  without  oscillating  wildly  between  gridpoints. 

We  can  state  these  goals  mathematically  as  follows.  We  seek  a  vector  Co  such  that 

7(co)  =  min/  (3.2) 

C 

where 

J(c)  =  ^  | n(rj)  -  n(rj)\2  +  a\  f  \h(r)  -  l(r)\2  dr  +  a2  f  | n'(r)  -  l'(r) |”  dr.  (3.3) 

j=  i  Jr0  Jiq 

Undefined  variables  are  given  in  the  following  table. 

We  do  not  put  a  weight  on  the  first  part  of  the  expression,  YiJ=i  \n(rj)  ~  n(rj)\  since  we  assume  that  the  weight  on  this 
part  of  the  expression  is  non-zero,  and  normalize  it  to  one.  Thus  a \  and  a2  are  decided  with  respect  to  the  importance 
of  the  first  term. 

Also,  we  choose  a\  >  a2,  since  the  second  and  third  terms  are  a  linear  combination  of  the  squared  L2  norm  and 
squared  Sobolev  norms.  That  is. 
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Table  3.1:  Table  of  Variables 


N 

the  number  of  points  at  which  we  have  data 

K 

the  dimension  of  the  spline;  the  number  of  basis  functions 

m 

the  order  of  the  spline 

iii 

the  smoothness  of  the  spline  at  the  knots 

n(rj) 

the  data  at  rj 

Kr ) 

linear  interpolation  between  the  known  gridpoints 

ai 

a  weight  that  determines  the  importance  of  our  approximation  h 
remaining  close  to  l(r) 

<*2 

a  weight  that  determines  the  importance  of  our  approximation’s 
derivative  h'  remaining  close  to  l'{r) 

O'!  r  | ti(r)  -  l(r) |2  dr  +  a2  f  | n{r)  -  l'(r) |”  dr 

Jr0  Jr0 


a\  || h{r)  -  l(r)\\2L2  +  a2  \\n{r)  -  l(r) ||2 


rrf  2  /  rrf 

a\  I  | h(r)  -  l(r)\  dr  +  ci2 

Jro 

rrf 

-l(r) |2  dr  +  | h\r)  -  l'(r)\ 

JrQ 

(fli  +  a2 )  f  | n(r)  -  I(r) |2  dr 

dr o 

+a2  I  |n'(r)  -  l'(r) |"  t/r. 

Jr0 


H 

f  I  n(r) 
Jr0 

A 


(3.4) 


Thus  we  have 


a  i  =  a\  +  a2  (3.5) 

=  A;  +  a2 

=>  a\  -  a2  =  ai  >  0.  (3.6) 

3.2  Necessary  and  sufficient  conditions  for  minimization 

The  necessary  condition  for  J  being  minimized  with  respect  to  c  is  that 

=0forA;=  (3.7) 

ock 

That  is,  for  k  =  1 . . .  K , 
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r*ff 

^  (n(rj)  -  n(rj) )  <t>k(jj)  +  at  I  (n(r)  -  /(r))  0*(r)  dr 
j=  i 

+a2  I  (n(r)-l'(r))(f>'k(r)dr 
Jro 


=  0 


N  f  K 

ZZ  Ci0i(O')  _  n(0') 

7=1  i=l 


ro 
7  X 


0*(r;)  +  a; 


Cj(pj(r)  -  l(r) 


+a2 


fZ 

v7r°  (=i 

rr/  K 

I  2  Q0j(r)  -  l'(r) 

Jro  ;=| 


0*(r)  Jr 


0*(r)  dr  =  0 


<f>i(r)<Pk(r )  Jr 


^  Q  ^  <pi(r j)<pk(r j)  -  n(rj)(pk(rj)  +  a,  ^  a  I 
i=l  V  7=1  J  y=l  1=1  ^r« 

ai  f  l(r)cpk(r)  dr  +  a2^Ci  f  0'(r)<^.(r)  Jr  -  a2  f  l'(r)<pk(r)  dr 
Jro  ;=1  Jro  Jr0 

f  rr/  r"  ,  , 

2^  <Pi(rj)<Pk(rj)  +  i  cf>i(r)(pk(r)  dr  +  a2  0;(r)0,(.(r)  dr 
/=1  l2=l  Jr»  Jr» 

\  ’  rrf  rrf 

-  Zj  n(rj)<pk(.rj)  -  a- 1  l{r)<pk{r)  dr  -  a2  l  (r)tf>k(r)  dr 


=  0 


=  0. 


Now,  expanding  this  for  all  A-  into  matrix  form. 


where 


(Ai  +  ffiA2  +  Q'2A3)c  =  (bi  +  aib2  +  a2b3) 


2li  MrjWiirj)  YlLi^iir^tpyiTj) 

i  MnWiirj)  I"=  i  <hirj»hirj ) 


’ZlLi<pK(rj)<pi(rj)  ' 
i  tKirjMrj) 


V  Ef=i  ‘t>i(rj)<pK(rj ) 

'  Xr/  0i(r)0i(r)  Jr 

£'  0 i(r)02(r)  dr 

,  £0f  <f>i(r)<pK(r)  dr 

'  XX  fiiirWiir)  dr 
XX  <p[(r)(f>'2(r)  dr 

,  XX^'i(r)^(r)c/r 


JP  020001  (r)  Jr 
fof  02(r)02(r)  Jr 


XX  02W0'i^)  ^ 

XX  02(r)02(r)  Jr 


EX=i  MrjWArj)  J 

£' f  </>K(r)<f>i(r)  dr  ' 
jf  0jsr(r)02(r)  Jr 

XX  <pK(r)(f>K(r)  dr  , 

XX  0^(r)0'i(r)  dr  ' 
XrX  PkWiW  dr 

XX  0i(r)0jrW  dr  > 


EX=1»(O)01(O)  1 

Ey=i  w(r y)02(r _,■) 

XX  l(r)<p\(r)dr  \ 

bi  = 

b2  = 

XX  ^(r)02(r)Jr 

,  Z/Li  n(rj)(/>K(rj)  , 

,  XX  l(r)<pk(r)dr  , 

(3.8) 


(3.9) 


(3.10) 


(3.11) 


(3.12) 


(3.13) 
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b3 


C  l'(r)4>\(r)dr  ' 
I/  l'{r)<p'2{r)dr 

C  l'(r)4>'k(j)dr  , 


Note  that  all  of  these  matrices  are  known  or  can  be  calculated.  Also  note  that  the  top  three  are  symmetric,  and  are 
generally  banded  since  B-splines  are  local  in  nature. 

The  sufficient  condition  is  that  J  must  be  convex  (concave  up)  with  respect  to  c.  We  know,  then,  that  the  Hessian  must 
be  positive  semi-definite.  The  Hessian  is  symmetric,  and  given  by 


H  =  Ai  +  A3  +  A3. 


(3.14) 


Hence  the  sufficient  condition  is  satisfied  if  the  eigenvalues  of  H  are  non-negative  [6],  The  matrices  Ai,  A3,  and  A3 
are  banded,  as  mentioned  above.  If  we  can  make  the  stronger  assumption  that  H  is  diagonally  dominant,  which  is 
probably  reasonable  for  B-splines,  the  Gershgorin  Circle  Theorem  [7]  states  that  all  real  parts  of  the  eigenvalues  will 
be  positive.  Further,  since  these  matrices  are  symmetric  (Hermitian)  H  is  also  Hermitian;  the  eigenvalues  are  therefore 
real. 

Hence,  as  long  as  the  sum  of  these  matrices  is  diagonally  dominant,  the  eigenvalues  will  be  positive,  and  the  sufficient 
condition  is  satisfied. 
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Chapter  4 

Approximation  in  two  dimensions 

4.1  The  objective  function  in  two  dimensions 

Recall  that  we  would  like  to  approximate  n(r,z),  a  function  on  two  variables  when  given  data  on  a  two-dimensional 
grid.  Letting  <p  represent  spline  basis  functions  in  the  r  direction,  and  <//  represent  spline  basis  functions  in  the  z 
direction,  we  seek  coefficients  ci;-  such  that 


K  P 

h{r,z )  =  EE  Cij<pi(r)if/j(z ) 
<= i  j= i 


(4.1) 


where  n  denotes  our  approximation  of  n,K  is  the  size  of  the  spline  basis  in  r  and  P  the  size  in  z.  Let  C  be  defined  by 


Cll 

C\2  ■' 

C\K 

C21 

C22  '  ' 

C2K 

-  Cp\ 

CPK 

Thus  our  spline  approximation  of  the  function  can  be  rewritten  in  matrix  form: 


(4.2) 


h(r,z)  =  xVT  C<I>  (4.3) 

where  'T  =  (t/q  (z), if p(z))T  and  ®  =  (cpi(r), (/>K(r))r. 

We  run  into  obstacles  when  we  move  to  two  dimensions  from  one,  however,  regarding  what  to  use  for  an  objective 
function.  Of  course  we  would  like  to  minimize  the  distance  between  our  approximation  and  our  data  at  points  where 
we  have  data;  this  is  not  hard  to  generalize.  However,  it  is  difficult  to  generalize  the  method  of  keeping  the  interpolating 
curve  smooth  between  grid  points. 

In  one  dimension,  the  linear  interpolation  between  grid  points  is  uniquely  defined  by  simply  connecting  the  data 
points  with  line  segments.  The  analog  to  this  in  two  dimensions  would  be  to  define  a  plane  between  sets  of  three 
points.  However,  we  are  working  with  a  rectangular  grid,  and  there  are  many  ways  to  choose  sets  of  three  points  that 
will  alter  our  approximation. 

In  Figure  4.1,  we  see  a  possible  problem  with  triangulating  the  grid  in  an  arbitrary  way.  If  we  divide  the  rectangles 
formed  by  the  grid  into  half  along  the  diagonal  we  see  that  one  interpolation  may  give  us  a  very  different  approximation 
for  points  not  on  the  grid  than  the  other,  though  this  example  is  extreme. 

In  Figure  4.1  we  see  two  different  triangulations  of  the  same  grid.  Call  the  first  t\  and  the  second  h_.  Using  a  lin¬ 
ear  combination  of  i\  and  U  is  one  possible  strategy  to  keep  our  approximation  from  oscillating  too  much  without 
restricting  it  in  an  artificial  manner.  We  define  the  objective  function  J  as  follows. 
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Figure  4.1:  Two  triangulations  of  a  rectangular  grid.  The  four  corner  points  are  hypothetical  data,  and  the  middle  point 
on  each  is  the  interpolation  in  each  case. 


m  =  X  X  \n(rj,Zj)  -  n(r;-,z7)|~  +  ai  I  I  |n(r,z)  -  ?i(r,z)|2  dz  dr 

i=  i  y=i  w^o  -'zo 

+  f  f  |n(r,z)  -  f2(r,z)|2  rfedrj  +  or2  (  f  f  \nz-t\J^  dz  dr  (4.4) 

Jz0  /  W'o  Jzo 

r'/  rz/,  ,2  r'/  rz/,  |2 

+  I  I  |«r  -  fi,r|  dzdr  +  I  |n,  -  f2,.-|  rfe  dr 

Jr0  Jzo  dr0  Jzo 

r'f  rzf ,  |2 

+  I  I  |«r  _  f2,r|  rfe  t/r 

Jr0  Jzo 

where  ?!  represents  the  planar  interpolation  where  the  grid  has  been  triangulated  in  one  direction  (see  Figure  4.1  for 
illustration)  and  f2  represents  the  other.  Note  that  the  definition  of  the  Sobolev  norm  tells  us  that  ct\  >  a2,  as  with  the 
1-D  case. 

Let  6,  y  be  integers  such  that  l<5<lVi,l<y<lV2  where  N\  is  the  number  of  grid  points  in  the  r  direction  and  /V2 
the  number  of  grid  points  in  the  z  direction.  Then 
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dJ  ^  ^2  Qfj,  /  f  r*^ f  dpi 

T —  =  2  V  V  (n(ritZj)  -  n(rhZj ))  - —  +  2  aJ  ( h(r,z )  -  h(r,z ))  7 —  dz  dr 

dcSy  jztjzt  1  '  dc6y  \JroJZo  dcSy 


r'rf('zf  dh 

I  (Mr,  z)  -  t2(r,  z))  77—  dzdr  |  +  2a2 

JroJzo 


dc , 


8y 


nZf 

<■ 

1 


h,  -  t\ ,)  -f—  c/s  c/r 
dc, 


<Sy 


ffflr  rnrzf  ^ 

(nr  -  fi,r)  7—^  dz  dr  +  (hz  -  t2tZ)  7—  c/s  dr 

Jc8y  JrnJzo  Jc8y 


=  2 


in 

Jr0Jzo  / 

a:  P  A^i  V2  A^i  V2 

zz  Ckp  Mn)<f>k(ri)  ^  ( ifrp(Zj)i//y(Zj ))  -  ^  0a(r/)  ^  n(rhZj)i//y(Zj) 

k=  1  P=1  /=1  7=1  i'=l  7=1 

f  r*rf  /'Z/  nzf 

2_J  2_J  CkP  I  <Pk(r)(ps(r)  I  (i//p(z)ifry(z))  dz  dr  -  I  <^(r)  I  (fi(r,s) 
A;=l  p=l  ^r0  Jzo  Jro  Jzo 

^  r*ry  r*Zf 

+h(r,  s))  iAy(s)  c/s  c/r)  +  a2  hZZ  Ckp  I  Mc)Mr)  I  (lAp(s)iAy(s))  c/sc/r 

,  A=1  p=i 

/r'A‘/  ^Z/  ^  ^  /r'r/ 

-  I  <M>')  I  (fi ,<:  +  f2,J  lAyfe)  c/s  c/r  +  2  ^  ^  qp  I  <p'k(r)<p’s(r) 

Jr0  Jzo  £=1  /;=1  Jr0 

f  (>t'p(z)i//y(zj)  dzdr-  f  0^(r)  f  (fLr  +  t2p)  if/y(z.)  dz  dr) 

Jzo  Jra  Jzn  / 


JZo  J r0  kJzo 

Our  necessary  condition  is  that  all  partials  are  zero.  That  is,  after  some  regrouping,  for  every  y,  <5, 


K  P  Ni  N2  K  P  r, 

zz  °kp  'Yj  Mn)Mn)  ^  MzjWvizj)  +  2<*i  ZZ  Ckp  I 

k=  1  p=  1  i=l  7=1  fc=l  /?=1  ^r° 


Mr)Mr) 


f 

Jzo 

K  P 


7= 

(  K  P 


il/p(z)i//y(z)  dz  dr  +  2a2 


iV  i  /^r/  /"*Z/ 

zz  Qp  I  <l>k(r)</>s(r)  I  'fk'pizM'y 

\k=l  p=  l  '''» 

\  AT]  M 


(s)  dzdr 


"  '  /^r/  /'Z/  "i 

ZZ  Ckp  I  (f>'k(r)(l>'s(r)  I  tl/p(z)if/y(.z)dzdr  -  ^jfo(ri)'^n(ri,Zj)i//y(Zj) 

k=  1  p=l  ^Zo  J  i=l  7=1 

-«I  f  f  T(r,zMy(z)  dz  dr  -  a2l  f  0'(r)  f  T  r(r,  z)if/y(z)  dz,  dr 

Jr0  Jzo  \Jro  Jzo 


rr/  r 

Mr) 

Jrn  Jzo 


Tz(r,z)i//'(z)  dz  dr  =  0 


(4.5) 


(4.6) 


where  T(r,z)  =  h(r,z )  +  t2(r,z). 

In  matrix  notation,  we  have  a  system  that  looks  like 

AiCDi  +  2aiA2CD2  +  2cr2  (A3CD2  +  A2CD3)  =  bi  +  crib2  +  cr2(b3  +  b4) 

where 


Ai 


MzjWliZi) 
2^i  MzpMzj) 


Z%Mzj)Mzj)  ' 

2^1  ilsp(Zj)i/kp(Zj)  , 


(4.7) 


(4.8) 
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■■■ 

EX‘i^i  (n)<f>K(n) 

Dr  =  : 

(4.9) 

,  E^i  Mri)<t>i(ri)  ■■■ 

EX‘i  MrdMo)  > 

XX  Ai(z¥t(z)  dz  ■■■ 

XX<Mz)(Mz)<fe  ' 

A2  =  | 

(4.10) 

,  XX  l/spizWlCz)  dz  ■■■ 

XX  <Mz)<Mz)  dz  , 

f^<t>i(r)tpi(r)  dr  ■■■ 

XX  <p\(r)<pK(r)  dr 

D2  =  : 

(4.11) 

,  XI'  <pK(r)<pi(r)  dr  ■■■ 

XX  <pK(r)<pK(r)  dr 

'  CKWUdz  ■■■ 

XX  *Ar  (z)«A^(z)  dz  ' 

A3  =  • 

(4.12) 

,  XX  V^(z)^(z)  dz  ■■■ 

XX  ^(z)'Ap(z)  *  , 

XX  (t>\  (r)4>\  (r)  dr  ■■■ 

XX  ^  'MMr)dr 

d3  = 

(4.13) 

.  .Q  <P'K(rWr)  dr  ■■■ 

XX  K^Mr)  dr 

lAlfel)  •••  Al(ZiV2)  \  ”^1,Zl-) 

n{rNl,Zi)  V  0l(ri 

)  •••  0/r(n)  ' 

(4.14) 

Mzi)  •••  Mzn2)  n(ruzN2)  : 

<Ov,,Za(2)  J'  ^1^1)  •••  ‘M'A)  J 

I  <P\(r)  3 

nr  1  T  nzf 

b2  =  I  :  T(r,z)(  <Ai(z),  •••,  Ap(z)  ) 

t/z  dr 

(4.15) 

Jr°  {  Mr)  r° 

n(  Mr)  )  r> 

ba  =  ;  Tz(  (z) 

•  •  •  ,  Ap(z)  )  t/z  c/r. 

(4.16) 

Ur0  Jz 0 

<Pk(j)  ; 

4>\  (r)  \ 

rrf  1  CZf  / 

b4  =  1  :  I  Tr[  (A, (z) 

•  •  •  ,  iAp(z)  )  &  c/r. 

(4.17) 

JrQ  Jzo 

(f>K(n  ; 

The  matrices  Ai,  D;  are  banded  (since  our  bases  are  B-splines)  and  symmetric  for  i  =  1, 2, 3;  A;,  I)j  as  well  as  b  can 
be  calculated.  However,  this  system  is  difficult  to  solve  in  its  current  form.  We  tackle  this  in  the  next  section. 

4.2  Solving  the  system 

We  replace  multiplication  of  three  matrices  with  multiplication  by  a  matrix  and  a  vector  to  make  the  system  easier  to 
solve.  That  is,  for  each  term  AjCDj, 

AiCDj  -»  A^c 

and  for  each  i 

bi  ->  b| 
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where  c  and  b'  are  vectors  of  length  k  x  p  and  A'„  is  a  square  matrix  with  k'/p  rows  and  columns.  Then  we  have  the 
linear  system 

(Aj  +  2a\A!2  +  2ai  (A^  A^J j  c  =  +  oqb)  +  ct 2  ^b^  +  b^j .  (4.18) 

The  modified  system  could  be  formed  in  different  ways,  but  for  consistency  we  describe  one  approach.  We  obtain  c 
by  stacking  up  the  columns  of  C,  one  on  top  of  the  other  from  left  to  right,  and  do  the  same  for  b[.  We  construct  the 
matrix  A'n  from  the  matrices  Ai  and  Dj  in  the  following  manner.  For  any  term  AjCDj  if 


then  we  let 


d\\  dn  •••  d\K 
dn  d22  ■  ■■  d2K 

l>i  -  ;  ;  ; 

dic\  .  diac 


(4.19) 


c/nAj 

<3^21  Ai 

^AlAi 

dnM 

d22  Ai 

tfcAi 

^lA'Aj 

^2A'Ai  •  • 

dKK^-i 

(4.20) 


Once  we  determine  the  values  for  c  from  Equation  4.18,  we  can  reassemble  c  back  into  its  original  matrix  form  C  by 
unstacking  the  columns,  and  use  Equation  4.3  to  find  an  approximation  of  the  function  n  at  any  point  (V,  z). 
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Chapter  5 

Tracing  rays  using  Runge-Kutta  methods 


We  can  apply  the  method  outlined  in  the  previous  chapters  to  determine  the  refractive  index  at  any  point  in  two 
dimensions.  Next  we  assume  that,  through  these  methods  or  analytically,  we  can  determine  the  refractive  index  at  any 
point  in  a  medium.  This  information  in  combination  with  the  ray  trace  equation  gives  us  an  approximation  for  the  path 
a  ray  may  travel  through  the  medium.  We  discuss  a  method  for  doing  this  in  this  chapter. 


5.1  The  ray  trace  equation 


Recall  that  the  ray  trace  equation  can  be  expressed  [5] 


d2r 

dt2 


—  nVn 


where  r  =  xi  +yj  +  zk,  and  n(r)  is  the  refractive  index  distribution  at  r.  t  is  defined  as 


(5.1) 


where  ds  is  an  arc  length  measure. 
Introducing  the  optical  ray  vector  T  as 


T  = 


dr 

dt 


allows  us  to  rewrite  the  second-order  ray  trace  equation  as  a  system  of  first-order  equations: 


(5.2) 


(5.3) 


dr 

dt 


=  T 


(5.4) 


—  =  nVn  (5.5) 

dt 

To  determine  the  path  of  a  ray  through  a  medium,  we  now  need  only  numerically  solve  this  system  of  ordinary 
differential  equations. 


5.2  Runge-Kutta  methods 

The  Runge-Kutta  methods  are  a  set  of  methods  that  numerically  approximate  the  solution  of  an  ordinary  differential 
equation  or  a  system  of  ordinary  differential  equations.  We  will  describe  the  idea  of  Runge-Kutta  methods,  the  classical 
Runge-Kutta  method,  and  finally  the  variation  we  chose  to  use  for  this  problem,  a  Runge-Kutta-Fehlberg  method 
known  as  RKF45. 
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5.2.1  General  Runge-Kutta  Methods 

Consider  the  initial  value  problem 


y'(t)  =  f(f,y),  yOo)  =  yo-  (5.6) 

Since  we  often  cannot  solve  such  a  system  analytically,  as  is  the  case  with  the  ray  equation,  we  use  the  information 
given  in  the  initial  value  problem  to  find  a  numerical  solution.  We  know  the  initial  value  of  the  vector  y,  and  we  know 
the  derivative  of  y  with  respect  to  t.  Using  the  derivative,  we  can  take  small  steps  along  t  starting  with  to  and  estimate 
y  at  each  step.  That  is,  for  a  reasonably  small  value  of  the  step  size  h, 

y(t  +  h)*y(t)  +  hf(t,y(t)).  (5.7) 

We  can  iteratively  use  the  above  expression  starting  with  the  initial  value  to  estimate  y  at  different  values  of  t.  This 
method  is  called  Euler’s  method  or  the  polygon  method  [6]. 

Euler’s  method,  however  is  a  method  of  order  1,  which  means  that  the  magnitude  of  the  difference  between  the  true 
difference  quotient  and  the  difference  quotient  for  the  approximation  is  0(hl)  =  0(h). 

Runge-Kutta  methods  are  one-step  methods  that  generalize  Euler’s  method  by  using  more  points  to  estimate  the 
function  value  at  each  step.  Using  a  Runge-Kutta  method  instead  of  Euler’s  method  allows  us  to  achieve  greater 
accuracy  with  larger  step  sizes.  For  example,  the  classical  Runge-Kutta  method  [6]  is  of  order  4.  The  classical  method 
is  given  by 


with 


y n  -  y„-i  +  7(ki  +  2k.2  +  2k.3  +  k4) 
o 


(5.8) 


kj  =  f  (/„_i,  y„_i) 


s  fl 

! 

h 

h  k, 

k i-i 

l  +  2,y"-1  +  ^r 

s  fl 

< 

h 

hk2 

k-i 

1  +  2,y"~1  +  T, 

k4  =  f(U_i  +  h,y  +  hkf). 


(5.9) 

(5.10) 

(5.11) 

(5.12) 


Note  the  similarity  to  Simpson’s  rule.  While  this  method  achieves  better  accuracy  than  Euler’s  method,  we  are  left 
with  the  problem  of  how  to  choose  h.  It  would  be  inefficient  to  use  trial  and  error  with  various  step  sizes,  trying  to 
determine  if  our  approximation  is  good  enough.  This  is  the  motivation  for  Runge-Kutta-Fehlberg  methods. 


5.2.2  Runge-Kutta-Fehlberg  methods 

Runge-Kutta-Fehlberg  methods  use  Runge-Kutta  methods  to  determine  whether  the  correct  step  size  h  is  being  used 
at  each  step,  and  to  choose  the  next  step  size  [6].  Specifically,  at  step  n  two  approximations  are  made:  one,  say  y„, 
using  a  Runge-Kutta  method  of  order  p,  and  the  other,  say  z„  a  Runge-Kutta  method  of  order  p  +  1 .  If  the  difference 
\yn  -  z„  \  is  below  a  certain  tolerance,  one  of  these  approximations  is  accepted,  a  step  size  for  the  next  step  is  calculated, 
and  the  procedure  is  repeated  for  the  next  step. 

Since  at  any  step  n  we  calculate  two  approximations,  y„  of  order  p  and  z.n  of  order  p  +  1,  we  much  choose  which  one 
to  use.  While  it  seems  logical  to  take  the  higher-order  approximation,  and  this  is  often  done,  the  error  analysis  done 
automatically  as  we  perform  the  Runge-Kutta-Fehlberg  method  applies  to  the  order  p  approximation.  It  is  therefore 
advisable  to  take  the  order  p  approximation  yn  particularly  in  the  case  of  stiff  problems.  [2] 
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5.2.3  Implementation 

We  employ  a  Runge-Kutta  pair  consisting  of  methods  of  order  4  and  5  as  given  in  [3].  To  move  from  y„  |  to  y„  we 
must  compute  the  following  six  vectors. 


ki 

k2 

k3 

k4 

k5 

kft 


1  (t/i-i ,  y„_i) 

f  [tn-  i  +  +  ^ki/z 


f  \tn-\  +  „h,yn-i  +  I  ^ki  +  — k2 1 h 


12 

t  ( tn-\  +  —  h,y„~i  + 


( 1932  7200,  7296, 

- kr - ki  + - k3  h 

\2197  2197  “  2197 


845 


,,  ,  439  0,  3680, 

I  t r  i  +  /z,  yn_i  +  - ki  —  8ki  + - k3 

1  J  1  \ 2 16  -  513  4104 


k Ah 


1  8  3544 

flU_1  +  -/z,y„_1+  --k1+2k2-— 1 


1859. 


11. 


-k4 - ks  /z 

4104  40 


(5.13) 

(5.14) 

(5.15) 

(5.16) 

(5.17) 

(5.18) 


Using  these  six  vectors,  two  approximations  are  made.  The  first  is  of  order  4: 


and  the  second,  of  order  5: 


y„  =  y„-i  + 


25  ,  1408, 

216kl  +  2565k3 


2197 

4104 


rksl/z 


(5.19) 


y«-i  + 


16 

135' 


ki  + 


6656,  28561, 

- k}  + - k4 

12825  56430 


h. 


(5.20) 


To  determine  whether  we  should  accept  one  of  these  approximations  for  the  mil  step,  we  test  whether  the  difference 
in  the  two  approximations  is  less  than  a  predetermined  error  control  tolerance,  e.  That  is,  we  accept  the  5th-order 
approximation  if 


|y„  -  z„|  <  e. 

The  value  of  h  for  the  next  step  is  then  chosen  by  finding  a  scalar  q  using  the  following  expression 

eh  \1/4 
^  I  Zn  —  y«  I  / 

We  determine  our  next  step  size,  say  hnew,  by  multiplying  q  with  h.  That  is 

knew  (5.23) 

Naturally  we  must  consider  the  possibility  of  the  denominator  of  q  being  0  when  y„  =  z„-  For  numerical  purposes, 
we  choose  a  maximum  step  size  value,  say  hmax,  so  that  if  at  any  time  q  is  very  large,  or  if  the  denominator  of  q  is  0, 
we  take  hmax  as  our  step  size  for  the  next  iteration.  We  also  choose  a  minimum  step  size,  /z„„„,  to  prevent  the  program 
from  becoming  too  expensive  to  run. 

5.2.4  Stability 

While  knowing  the  order  of  the  Runge-Kutta  method  we  use  gives  us  an  estimate  in  terms  of  h  of  the  order  of  the 
error  per  step  in  our  method,  we  still  must  keep  in  mind  the  possibility  of  instability  in  our  method,  leading  to  an 
approximation  of  the  solution  to  the  differential  equation  that  grows  further  away  from  the  solution  as  t  increases. 

To  consider  the  possibility  of  stability  problems,  we  find  our  region  of  absolute  stability  in  the  complex  plane.  Consider 
the  ordinary  differential  equation 


(5.21) 


(5.22) 
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Figure  5.1:  Although  the  solution  approaches  0  as  t  — »  oo,  the  Euler  method  approximation  using  step  size  h  -  1/7 
does  not.  This  is  an  example  of  instability. 


y'(t)  =  Ay(t).  (5.24) 

We  can  express  a  Runge-Kutta  method  in  vector  form  as  follows  [1], 

Yi  =  y„-\  (5.25) 

Y2  =  y„-i  +  ha2\f(Yx)  (5.26) 

Ys  =  yn-\  +  h\as\f{Y{)  +  aS2f{Y2)  +  ■  ■  ■  +  aS'S-\f(Ys-\)],  (5.27) 

y„  =  yn-i+h[blf{Yl)  +  b2f{Y2)  +  ---b2f{Ys)\  (5.28) 


In  the  fifth-order  method  described  in  (5.20),  for  example,  k;  =  f(Y ,)  for  i  =  1, 2, . . . ,  s  =  6.  This  can  equivalently  be 
written: 


Y=y„_,e  +  /iAf(Y). 

where  Y  =  [Y\,Y2, YS]T ,  e  =  [1, 1, . . . ,  l]r  and 


A  = 


0  ■■■ 

a2\  0 

«3i  a22 


0 


,  1 


...  O' 

...  o 
...  o 

0 


(5.29) 


(5.30) 
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Also,  let  b  =  [b\,b2,.  bs]T  and  let  z  =  Ah.  Using  the  properties  of  the  simple  ODE  (1)  along  with  (5.29)  gives  us 


Y  =  y„_je+zAY 
y„  =  yn-\  +  zbrY 

We  would  like  to  find  the  region  of  stability.  From  [  1  ]  the  function  r(z)  determining  this  is 


(5.31) 

(5.32) 


yn 

y„- 1 


(5.33) 


=  1  + 


zb'Y 

yn-\ 


(using  Equation  5.32) 


=  1  +  zb 


yn-i 


=  1  +  "br  (e  +  — AY)  (using  Equation  5.31) 

>’«- 1 


=  1  +  zb' 


(  0  . 

«21  0 

031  O32  0 


y„- 1 


••  0  1 

0 

■■  0 
•.  0 

V  osi  .  oStS~ i  0 

=  1  +  zbr  (i +zA +z2A2  +  ■  •  • +z^'Ai_1)e 

For  a  Runge-Kutta  method  of  order  p,  if  k  <  p 


yn- 1 

y,,~  i  +  ha2  \f(Y\) 

}’n- 1  +  h  \as\f(Yi))  +  •  •  •  +  fls,s-i/(Es-i)] 


from  [1].  Thus 


r(z)  =  l+  z+  2-  +  --  -  +  :-  +  cp+izp+l  +  •  •  •  +  cszs 


2! 


P- 


(5.34) 


(5.35) 


(5.36) 


where  for  i  =  p  +  l,  p  +  2, ...  s,  the  coefficient  c,-  =  Ir  A'  'e.  Now,  for  our  situation  we  have 


0 

0 

0 

0 

1 

4 

0 

0 

0 

3 

32 

9 

32 

0 

0 

1932 

7200 

7296 

0 

2197 

2197 

2197 

439 

-8 

2 

3680 

845 

216 

^44 

ii^4 

~  27 

2565 

4104 

0  0  ' 
0  0 
0  0 
0  0 
0  0 


and 


(5.37) 


We  therefore  have  the  polynomial 


br 


16  6656  28561  -9  2 

135 ’ 12825’  56430'  50’  55 


(5.38) 


,2  73  ,4  5 

r(z)  =  1+  z+  Y  +  ^  +  ^  +  Y^  +  -00048076923077z6.  (5.39) 

To  determine  the  stability  region,  we  must  find  out  when  r(z)  <  1.  Using  a  routine  from  [1]  we  plotted  the  stability 
region  in  the  complex  plane,  the  interior  of  which  is  the  set  of  values  of  z  for  which  the  approximation  of  the  test 
function  y'(t)  -  Ay(t )  is  stable. 
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Stability  Region  for  RK4  and  RK5 

4 - i - i - i - - 


-41 - * - ‘ - 1 - - 

-4  -3  -2-1  0  1 

Re(z) 


Figure  5.2:  The  region  of  absolute  stability  for  RK4  and  RK5  given  in  (5.19)  and  (5.20). 


While  we  use  a  Runge-Kutta-Fehlberg  method  rather  than  a  Runge-Kutta  method,  our  approximations  are  in  fact 
determined  by  the  Runge-Kutta  method  of  order  4  or  5  whose  stability  region  is  shown.  Therefore,  assuming  li  is 
always  chosen  such  that  z  =  All  is  in  the  intersection  of  the  two  stability  regions,  our  method  is  also  stable. 
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Chapter  6 

Conclusion 


Given  data  values  for  a  function  on  a  two-dimensional  rectangular  grid,  using  B-splines  allows  us  to  estimate  with 
relative  computational  ease  reasonable  values  for  the  function  at  any  point  in  the  space  covered  by  the  grid.  Two 
advantages  that  B-splines  provide  are  their  local  nature,  and  the  ability  to  combine  them  using  tensor  products.  We 
demonstrated  in  this  report  how  to  make  the  system  of  tensors  that  results  from  using  a  B-spline  basis  in  two  dimensions 
into  a  more  standard  matrix  system. 

The  next  step  is  to  extend  the  use  of  a  B-spline  basis  into  three  dimensions,  and  then  to  n  dimensions  for  any  integer 
n  >  1.  For  a  cylindrically  symmetric  three-dimensional  space,  the  two-dimensional  procedure  that  we  have  outlined 
should  be  sufficient  if  the  problem  is  set  up  to  exploit  the  symmetry.  For  a  more  general  space,  extension  to  n  uses  the 
same  idea  as  the  problem  in  two  dimensions,  and  in  theory  is  not  much  more  difficult.  The  challenge,  however,  lies  in 
restructuring  a  system  of  high-rank  tensors  into  a  system  which  can  be  implemented  more  easily  on  a  computer. 
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Appendix  A 

Code 


Given  data  on  a  two-dimensional  grid,  we  implement  the  theory  mentioned  in  the  previous  sections  to  approximate 
the  function  of  interest.  We  develop  the  following  routines  to  achieve  this.  An  implementation  in  C  is  given  in  A.l, 
followed  by  an  example.  MATLAB  code  is  given  in  A. 3. 


A.l  C  code 

A.l.l  VBASIS 

#include  "spline. h" 

output  vBasis(double  x[] ,  double  grid[] ,  int  order,  int  derivative,  int  xLength, 
int  gridLength) 

{ 

*  ©Author:  Dane  Burrows 

*  ©Date  9-July-07 

*  ©Description: 

*  This  routine  evaluates  the  values  of  the  B-spline  basis 

*  (or  the  derivatives)  functions  at  given  points.  The  grid 

*  points  and  the  order  of  the  splines  are  specified  by  the 

*  user,  however,  additional  grid  points  outside  of  the 

*  interval  [xmin,  xmax]  are  chosen  by  the  program  to  provide 

*  a  complete  basis. 

*  ©Usage: 

*  output  <name>=vBasis(x,  grid,  order,  derivative,  xLength,  gridLength); 

*  Input : 

*  x  :  array  of  values  for  x  on  which  the  basis  functions 

*  are  to  be  evaluated. 

*  grid  :  the  grid  points  in  ascending  order,  all  grid  points  must 

*  be  distinct.  The  interval  on  which  the  spline  basis  functions 

*  are  defined  are  given  by: 

*  [grid[®] ,  grid[N] ]  . 
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*  where  N  is  the  length  of  the  array  grid. 

*  order  :  order  of  the  spline  functions. 

*  derivative  :  order  of  derivative  needed. 

*  xLength  :  an  integer  value  showing  the  length  of  the  array  x. 

*  gridLength  :  an  integer  value  showing  the  length  of  the  array  grid. 

*  Output : 

*  v  :  an  array  of  dimension  order  +1  by  M,  where  M  is  the  length  of 

*  the  array  x. 

*  ndim  :  total  number  of  basis  elements,  ndim=N+order- 1 . 

*  index  :  indices  of  the  basis  elements  with  non-zero  values  at  a 

*  point  x.  index  is  a  2  by  M  array, 

*  index [®] [k]  --  lowest  index  of  non-zero  basis  element  at  x[k] . 

*  index [1] [k]  --  highest  index  of  non-zero  basis  element  at  x[k] . 

*  ©Note: 

*  Output  is  a  structure  defined  in  functions. h. 


output  out; 

int  i,  j,  k,  factor,  lcount=®,  rcount=®,  acount=®,  counter=®,  M=xLength, 
N=gridLength,  lgridLength=order ,  rgridLength=order ; 

double  localg [2*order+l] ,  trunc [2*order+2] ,  dgrid,  n,  lgrid[lgridLength] , 
rgrid[rgridLength] ,  agrid[lgridLength+gridLength+rgridLength] ; 

out. ndim  =  N+order-1; 

out . index=Array2D<int> (2 ,  M) ; 

out.v=Array2D<double>(order+l,  M) ; 

//Construct  the  augmented  grid 

// 

dgrid=grid[l] -grid[Q] ; 
for (i=® ; i<lgridLength ; i++) 

{ 

lgrid[i]=dgrid*i  +  grid[®] -order*dgrid; 
lcount++ ; 

} 

dgrid=grid[N-l] -grid[N-2] ; 
for (i=® ; i<rgridLength ; i++) 

{ 

rgrid[i]=grid[N-l]  +dgrid*(i+l) ; 
rcount++ ; 

} 

for (i=® ; i<lcount ; i++) 

{ 

agrid[acount]=lgrid[i]  ; 
acount++ ; 

} 

for (i=® ; i<gridLength; i++) 

{ 

agrid[acount]=grid[i] ; 
acount++ ; 
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} 

for (i=0 ; i<rcount ; i++) 

{ 

agrid[acount]=rgrid[i] ; 
acount++ ; 


//Main  loop  over  points  x 

// 

for(k=0;  k<M;  k++) 

{ 

for(j=0;  j<N-l;  j++) 

{ 

if((sign(x[k]-grid[j])*sign(grid[j+l]-x[k]))>=0) 

break; 

} 

if (x[k]<grid[0] ) 

{ 

j=«; 

} 

if(x[k]>grid[N-l]) 

{ 

j=N-2 ; 

} 

//Evaluate  the  values  of  the  basis  functions  (or  derivatives)  at  x(k) 

// 

//  1.  Evaluate  the  values  of  the  truncated  polynomials 

// 

factor=l ; 
if (derivative  >0) 

{ 

for(i=0;  i<=derivative-l ; i++) 

{ 

factor=factor*(order-i) ; 

} 

} 

for (i=0 ; i<2*order+2 ; i++) 

{ 

trunc[i]=0; 

} 

counter=0 ; 

for(i=j ;i<=j+2*order+2;i++) 

{ 

localg [counter]=agrid[i]  ; 
counter++; 

} 

for (i=0 ; i<order*2+2 ; i++) 

{ 

if(0<=(x[k] -agrid[j+i] )) 
trunc[i]=x[k]-agrid[j+i] ; 
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if(order  >  derivative) 

trunc [i]=factor*(pow(trunc[i] ,  order-derivative)) ; 
else  if(order  ==  derivative) 
trunc [i]=factor*sign(trunc [i] ) ; 
else 
{ 

cout  «  "The  spline  function  is  not  differentiable"; 
exit(-l) ; 

} 

} 

111.  Compute  the  divided  differences 

// 

int  1,  11; 

for (1=0 ; l<order+l ; 1++) 

{ 

for (11=0;  ll<2*order+2-l ;  11++) 

{ 

double  tmp=(trunc [11+1] -trunc [11] )/(localg [11+1+1] -localg [11] )  ; 
trunc [11] =tmp; 

} 

} 

//3 .  Store  the  value  in  the  vector  v 

// 

double  itrunc[2] [1] ; 
itrunc [0] [0]=j ; 
out . index [0] [k] = j ; 
itrunc[l] [0]=j+order; 
out. index [1] [k]=j+order; 
for (i=0 ; i<order+l ; i++) 

{ 

out.v[i] [k]=trunc[i] ; 

} 


return  out; 

} 

A.1.2  VBNEUMANN 

#include  "spline. h" 

output  vbneumann(double  x[],  double  grid[] ,  int  order,  int  derivative, 
int  xLength,  int  gridLength) 

{ 


*  ©Author:  Dane  Burrows 

*  ©Date  9-July-07 
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*  ©Description: 

*  Evaluuates  the  value  of  the  basis  elements  of  spline  functions  of  the 

*  specified  order  on  the  given  grid  which  satisfies  the  Neumann  boundary 

*  conditions. 

*  ©Usage: 

*  output  <name>=vbneumann(x,  grid,  order,  derivative,  xLength, 


gridLength) ; 

Input : 

X 

:  array  of  values  for  x  on  which  the  basis  functions 
are  to  be  evaluated. 

grid 

:  the  grid  points  in  ascending  order,  all  grid  points 
must  be  distinct.  The  interval  on  which  the  spline 
basis  functions  are  defined  are  given  by: 

[grid[0] ,  grid[N]]. 

where  N  is  the  length  of  the  array  grid. 

order 

:  order  of  the  spline  functions. 

derivative 

:  order  of  derivative  needed. 

xLength 

:  an  integer  value  showing  the  length  of  the  array  x. 

gridLength 

:  an  integer  value  showing  the  length  of  the  array  grid 

Output : 

V 

:  an  array  of  dimension  order  +1  by  M,  where  M  is  the 
length  of 
the  array  x. 

ndim 

:  total  number  of  basis  elements,  ndim=N+order- 1 . 

index 

:  indices  of  the  basis  elements  with  non-zero  values  at 
point  x.  index  is  a  2  by  M  array, 

index[0][k]  --  lowest  index  of  non-zero  basis 
element  at  x[k] . 

index[l][k]  --  highest  index  of  non-zero  basis 

*  element  at  x[k] . 

*  ©Note: 

*  Output  is  a  structure  defined  in  functions. h. 


Array2D<double>  tau,  v,  u; 

Array2D<int>  index; 
output  tmp,  out; 

int  M=xLength,  N=gridLength,  k,  i,  ndim; 
double  interval [2]; 

interval [0]=grid[0] ; 
interval [l]=grid[N-l] ; 

tmp  =  vBasi s (interval ,  grid,  order,  1,  2,  gridLength); 
tau=tmp . v ; 

tmp  =  vBasis(x,  grid,  order,  derivative,  xLength,  gridLength); 

ndim=tmp.ndim; 

u=tmp . v ; 
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v=u. copyO  ; 
index=tmp . index ; 

for(k=0;  k<M;  k++) 

{ 

if(index[0] [k] !=1) 

index [0] [k]=index[0] [k] ; 
else 
{ 

for (i=0 ; i<order ; i++) 

v[i] [k]=u[i+l] [k]-tau[i+l] [0]*u[0] [k]/tau[0] [0] ; 

} 

if(index[l] [k] !=ndim) 

index [1] [k]=index[l] [k] -1 ; 
else 
{ 

for (1=0 ; i<order ; i++) 

v[i] [k]=v[i] [k]-tau[i] [1] *u[order+l] [k]/tau[order+l] [1] ; 
index [1] [k]=ndim-2; 

} 

ndim=ndim-2 ; 
out . ndim=ndim ; 
out . v=v ; 

out . index=index ; 
return  out; 

} 

} 

A.1.3  LEAST  SQUARES  APPROXIMATION 

#include  "spline. h" 

Array lD<double>  lsqapp(ArraylD<double>  xdata,  Array lD<double>  ydata, 
Array lD<double>  wdata,  ArraylD<double>  xgrid,  int  order,  double  alphaO, 
double  alphal) 

{ 


*  ©Author:  Dane  Burrows 

*  ©Date  9-July-07 

*  ©Description: 

*  Compute  the  least  square  approximation  of  the  data  set  using  a  given 

*  set  polynomial  spline  functions.  The  optimization  functional  is  given 

*  by 

*  J(coef)  =  \sum"N_{J=l}  wdata_j  |L(t_j)-S(t_j)  |  "2 

*  \alpha_0\int''{t_max}_{t_min}  I  L[t] -S[t]  |  ~2dt 

*  \alpha_l\int''{t_max}_{t_min}  | L ’  [t]-S’  [t]  |~2dt, 

*  where: 
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N  :  number  of  data  points. 

L  :  liner  spline  interpolation  of  the  data. 
S  :  polynomial  spline  function. 

\alpha_®  :  where  on  the  L_2  norm. 
\alpha_l  :  weight  on  the  H_1  norm. 


*  ©Usage: 

*  Array lD<double>  <name>=vBasis(xdata,  ydata,  wdata,  xgrid,  order, 

*  alpha®,  alphal) ; 


Input : 

xdata  :  data  values  for  the  independent  variable. 

ydata  :  data  values  for  the  dependent  variable. 

wdata  :  weights  on  the  data  points. 

xgrid  :  grid  for  the  spline  function. 

order  :  order  of  the  polynomial  spline. 

alpha®  :  weight  on  the  L_2  norm. 

alphal  :  weights  on  the  H_1  norm. 


Output : 

coef  :  coefficients  for  the  optimal  spline  function. 


output  start=vBasis (xdata,  xdata,  1,  ®,  xdata. dim() ,  xdata. dim()); 
Array2D<double>  PI (start. ndim, start. ndim) ,  P2 ,  P3,  W(wdata . dim() , 
wdata. dim()),  Al,  A2,  A3,  A,  Q,  intp_tmp; 

Array lD<double>  rl,  r2 ,  r3,  r,  coef,  intp; 
double  xmin,  xmax; 
int  i,  j,  k; 

//Evaluate  the  pointwise  term. 

// 

for (i=® ; i<start . ndim ; i++) 

{ 

for ( j  =® ; j  <start . ndim ; j  ++) 

{ 

Pi C j ] [i]=®; 

} 

} 

for(k=Q;  k<start .ndim;  k++) 

{ 

for(j=start . index[Q] [k] ;  j<=start . index[l] [k] ;  j++) 

{ 

P 1 [k] [ j ]  =  start. v[j-start. index[0] [k]] [k] ; 

} 

} 

intp=inverse(Pl)*ydata; 

output  filter=vBasis(xdata,  xgrid,  order,  ®,  xdata. dim() ,  xgrid. dim()); 
Q=Array2D<double>  (start. ndim,  filter .ndim) ; 
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for(int  k=® ;k<start .ndim;  k++) 

{ 

for(int  j=filter . index[®] [k] ;  j<=filter . index [1] [k] ; j++) 

{ 

Q[k] [j]=filter.v[j-filter.index[®] [k] ] [k] ; 

} 

} 

for(i=®;i<wdata.dim() ;i++) 

{ 

for(j=Q ; j<wdata. dim() ; j++) 

{ 

W[j]  [i] =® ; 

} 

} 

for(i=®;i<wdata.dim() ;i++) 

{ 

W[i] [i]=wdata[i] ; 

} 

for(i=Q;i<intp_tmp.dimlO ;i++) 

{ 

for(int  k=®;  k<intp_tmp . dim2 () ;  k++) 

{ 

intp[i]+=intp_tmp[i] [k]*ydata[k] ; 

} 

} 

Al=Array2D<double>  (filter .ndim,  filter .ndim) ; 
r l=transpose (Q) *W*ydata ; 

Al=transpose (Q) *W*Q ; 

//Evaluate  the  L_2  term 

// 

xmin=xdata[®] ; 
for(i=l;i<xdata.dim() ;i++) 

{ 

i f (xdata [i ] <xmin) 
xmin=xdata[i] ; 

} 

for(i=®;i<xgrid.dim() ;i++) 

{ 

i f (xdata [i ] <xmin) 
xmin=xgrid[i] ; 

} 

xmax=xdata[Q]  ; 
for(i=l;i<xdata.dim() ;i++) 

{ 

if(xdata[i] >xmax) 
xmax=xdata[i] ; 

} 

for(i=®;i<xgrid.dim() ;i++) 

{ 

if(xdata[i] >xmax) 
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xmax=xgrid[i]  ; 


} 

A=Array2D<double>  (filter .ndim,  filter .ndim) ; 

P2=Array2D<double>  (xgrid.dimQ,  xdata.dimO); 

A2=Array2D<double>  (filter .ndim,  filter .ndim) ; 

A3=Array2D<double>  (filter .ndim,  filter .ndim) ; 

P3=Array2D<double>  (xgrid.dim() ,  xdata.dimO); 

P2=innprd(xgrid,  order,  0,  xdata,  1,  0,  xmin,  xmax,  xgrid.dim(), 
xdata. dim()) ; 
r2=P2*intp ; 

A2=innprd(xgrid,  order,  0,  xgrid,  order,  0,  xmin,  xmax,  xgrid.dim(), 
xgrid. dim()) ; 

//Evaluate  the  H_1  term 

// 

P3=innprd(xgrid,  order,  1,  xdata,  1,  1,  xmin,  xmax,  xgrid. dim(),  xdata.dimO); 
r3=P3*intp ; 

A3=innprd(xgrid,  order,  1,  xgrid,  order,  1,  xmin,  xmax,  xgrid. dim(), 
xgrid. dim()) ; 

//Solve  for  the  optimal  coefficients 

// 

r=rl+alpha0*r2+alphal*r3 ; 

A=Al+alpha0*A2+alphal*A3 ; 
coef  =  inverse (A) *r; 

return  coef; 

} 

A.1.4  V  SPLINE 

#include  "spline. h" 

Array lD<double>  vspline (Array lD<double>  x,  ArraylD<double>  grid, 
int  order,  int  derivative,  Array lD<double>  coef) 

{ 


*  ©Author:  Dane  Burrows 

*  ©Date  9-July-07 

*  ©Description: 

*  Evaluate  a  given  polynomial  spline  function. 

*  ©Usage: 

*  Array2D<double>  <name>=vspline(x,  grid,  order,  derivative,  coef); 

*  Input : 

*  x  :  values  of  the  independant  variable. 

*  xgrid  :  grid  on  which  the  splines  are  defined. 
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*  order  :  order  of  spline. 

*  dev  :  order  of  derivative. 

*  coef  :  coefficients  with  respect  to  the  standard  basis. 

*  Output : 

*  v  :  value  of  the  spline. 


Array lD<double>  v(x.dimO); 
output  tmp; 
int  i,  j; 

tmp=vBasis(x,  grid,  order,  derivative,  x.dimO,  grid.dimQ); 

if (tmp .ndim! =coef . dim()) 

{ 

cout  «  "The  dimension  of  the  coefficient  vector  is  wrong.  " 

«  tmp. ndim  «  "  "  «  coef.dimO  «  endl; 
exit(-l)  ; 

} 

//Calculate  the  spline  values. 

// 

for(i=Q;  i<x.dim();  i++) 

{ 

v[i]=®; 

for(j=tmp . index[0] [i] ;  j<=tmp . index [1] [i] ;  j++) 

{ 

v[i]=v[i]+coef[j]*tmp.v[j-tmp.index[Q] [i] ] [i] ; 

} 

} 

return  v; 

} 

A.1.5  2D-SLICE  COEFFICIENT 

#include  "spline. h" 

Array2D<double>  slice_coef (Array lD<double>  z,  Array lD<double>  nr, 
Array2D<double>  r,  Array2D<double>  nval,  Array lD<double>  zgrid, 

Array lD<double>  nrgrid,  Array2D<double>  rgrid,  int  order,  double  alpha®, 
double  alphal) 

{ 


*  ©Author:  Dane  Burrows 

*  ©Date  9-July-®7 
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©Description: 

Evaluate  the  coefficient  for  each  slice  at  p_{ i }  to  approximate  the  data  set 
r(i,l:nr[i])  and  nval(i , 1 :nr [i] ) .  This  approximation  is  done  using  one 
dimensional  approximation. 

©Usage: 

Array2D<double>  <name>=slice_coef(z,  nr,  r,  nval,  zgrid,  nrgrid,  rgrid, 
order,  alpha®,  alphal); 

Input : 

zgrid  :  grid  points  in  z. 

nrgrid  :  number  of  grid  points  in  r  at  each  slice. 

rgrid  :  grid  points  in  r. 

order  :  order  of  the  spline  requested. 

z  :  z  values  for  data. 

nr  :  number  of  r  data  points  at  each  slice. 

r  :  r  measurements. 

nval  :  intensity  measurements. 

alpha®  :  weight  for  data  approximation. 

alphal  :  weight  for  derivative  approximation. 

Output : 

coefz  :  lenz  by  (nt[i]+order-l)  2  dimensional  array  containing 
the  approximation  coefficients  at  each  slice. 


int  lenz=z.dim(),  i,  j; 

Array2D<double>  coefz (lenz,  (int)nr[0]+order-l) ; 

for(i=®;  i<lenz;  i++) 

{ 

Array lD<double>  tmp,  tmpr,  tmpnval,  ones,  tmprgrid; 
tmpr=ArraylD<double>  ((int)nr[i]) ; 
tmpnval=ArraylD<double>  ((int)nr[i]) ; 
ones=ArraylD<double>  ((int)nr[i]) ; 
tmprgrid=ArraylD<double>  ((int)nrgrid[i]) ; 

for (j=® ;  j<nr[i] ;  j++) 

{ 

tmpr [j]=r [i] [j] ; 
tmpnval [j]=nval [i]  [j]  ; 
ones [ j ] = 1 ; 

} 

for(j=®;  j<nrgrid[i] ;  j++) 

{ 

tmprgrid [j ] =rgrid [i] [j] ; 

} 

tmp=lsqapp(tmpr ,  tmpnval,  ones,  tmprgrid,  order,  alpha®,  alphal); 
for(j=®;  j<tmp.dim();  j++) 
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{ 

coefz [i] [j]=tmp[j] ; 

} 


return  coefz; 

} 


A.1.6  INNER  PRODUCT 

#include  "spline. h" 

Array2D<double>  innprd(double  grdl[],  int  ordl,  int  devl,  double  grd2[], 
int  ord2,  int  dev2 ,  double  xmin,  double  xraax,  int  grdlLength,  int  grd2Length) 

{ 


*  ©Author:  Dane  Burrows 

*  ©Date  9-July-Q7 

*  ©Description: 

*  Computes  the  matrix  of  the  inner  product  of  two  families  of  polynomial 

*  spline  basis  functions.  If  {B_K} ,  k=l .  N  and  {C_j},  j=l,  ...,  M, 

*  then  the  matrix  is  given  by: 

*  A_{k,  j}=<B_k,  C_j>. 

*  ©Usage: 

*  Array2D  <name>=innprd(grdl ,  ordl,  devl,  grd2 ,  ord2 ,  dev2 ,  xmin,  xmax) ; 

*  Input : 

*  grdl  :  grid  of  points  of  the  first  group  of  polynomial  spline 

*  functions. 

*  ordl  :  the  order  of  the  first  group  of  spline  functions. 

*  devl  :  order  of  the  derivatives  of  the  first  group  of 

*  spline  functions. 

*  grd2  :  grid  points  of  the  second  group  of  polynomial  spline 

*  functions. 

*  ord2  :  the  order  of  the  second  group  of  spline  functions. 

*  dev2  :  order  of  the  derivatives  of  the  second  group  of 

*  spline  functions. 

*  xmin  :  lower  bound  of  the  interval  of  integration. 

*  xmax  :  upper  bound  of  the  interval  of  integration. 

*  Output : 

*  A  :  the  matrix  of  the  inner  product. 

*  ©Note: 

*  Uses  the  wt  function  to  determine  the  weights  (needs  to  be  changed 

*  if  values  other  than  1  are  desired) . 
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int  z=®,  nint=5,  i,  j,  k=®,  n=grdlLength,  lgridLength=ordl ,  acount=®, 
rgridLength=ordl ,  Nl,  N2,  icount=Q; 
double  weight[5],  x[5] ,  dgrd=grdl [1] -grdl [®] ,  lgrd[lgridLength] , 
rgrd[rgridLength] , 

agrdl [lgridLength+grdlLength+rgridLength] ,  gridl[l],  grid2[l],  u; 
double  alpha[5] ,  y [ 5] ; 

Array2D<double>  A; 
output  fltrl,  fltr2; 

weight [®]=® .2369268851 ; 
weight [1]=® . 47862 867® 5 ; 
weight[2]=®. 5688888889; 
weight[3]=weight[l] ; 
weight [4] =weight[®] ; 
x[®]=-®. 9061798459; 
x[l]=-®. 5384693101; 
x[2]=0; 
x[3]=-x[l] ; 
x[4]=-x[0] ; 

//construct  the  combined  grid 

// 

for (i=® ; i<lgridLength ; i++) 

{ 

lgrd[i]=dgrd*i  +  grdl [®] -ordl*dgrd; 

} 

dgrd=grdl[n-l]  -  grdl[n-2]; 
for (i=® ; i<rgridLength ; i++) 

{ 

rgrd[i]=grdl [n-1]  +  dgrd*(i+l) ; 

} 

for (i=® ; i<ordl ; i++) 

{ 

agrdl [acount] =lgrd[i] ; 
acount++ ; 

} 

for (i=® ; i<grdlLength; i++) 

{ 

agrdl [acount] =grdl[i] ; 
acount++ ; 

} 

for (i=® ; i<ordl ; i++) 

{ 

agrdl [acount] =rgrd[i] ; 
acount++ ; 

} 

lgridLength=ord2 ; 
rgridLength=ord2 ; 
n=grd2Length; 
acount=® ; 
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dgrd=grd2 [1] -grd2 [0] ; 
for (i=® ; i<lgridLength ; i++) 

{ 

lgrd[i]=dgrd*i  +  grd2 [®] -ord2*dgrd; 

} 

double  agrd2 [lgridLength+grd2Length+rgridLength] ; 
dgrd=grd2 [n- 1] -grd2 [n-2] ; 
for (i=® ; i<rgridLength ; i++) 

{ 

rgrd[i]=grd2 [n-1]  +dgrd*(i+l) ; 

} 

for (i=® ; i<ord2 ; i++) 

{ 

agrd2 [acount] =lgrd[i] ; 
acount++; 

} 

for (i=® ; i<grd2Length ; i++) 

{ 

agrd2 [acount]=grd2 [i] ; 
acount++; 

} 

for (i=® ; i<ord2 ; i++) 

{ 

agrd2 [acount]=rgrd[i]  ; 
acount++; 

} 

double  cgrid[ordl*2+ord2*2+grdlLength+grd2Length] ; 
int  ccount=®; 

for (i=® ; i<ordl* 2+grdlLength ; i++) 
cgrid[ccount++]=agrdl [i]  ; 
for (i=® ; i<ord2* 2+grd2Length ; i++) 
cgrid[ccount++]=agrd2 [i] ; 
int  elements=  sizeof (cgrid)/sizeof (double) ; 
sort(cgrid,  elements+cgrid) ; 

double  igrd[2+ccount] ; 
igrd [icount++] =xmin ; 
double  cx=cgrid[8] ; 

i-® ; 

while (cx !=cgrid[ccount-l] ) 

{ 

if(cx  >=  xmax) 
break; 

if (igrd[k]<cx) 

{ 

igrd [icount++] =cx ; 
k++; 

} 

cx=cgrid[i++] ; 

} 

igrd [icount++] =xmax ; 
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gridl [®]=grdl [®] ; 
grid2 [Q]=grd2 [®] ; 

fltrl=vBasis(gridl ,  grdl,  ordl,  devl,  1,  grdlLength) ; 
fltr2=vBasis(grid2 ,  grd2,  ord2 ,  dev2 ,  1,  grd2Length) ; 

//Calculate  the  inner  product  matrix 

// 

Nl=fltrl .ndim; 

N2=fltr2 .ndim; 

A=Array2D<double>(Nl ,N2)  ; 
for(i=® ; i<fltr2 .ndim; i++) 

{ 

for(j=Q ; j<fltrl .ndim; j++) 

{ 

A[j][i]=8; 

} 

} 

for (z=® ; z<icount- 1 ; z++) 

{ 

double  a=igrd[z] ; 
double  b=igrd[z+l] ; 
for(i=Q ; i<5 ; i++) 

y [i]=(b+a)/2+(b-a)*x[i]/2 ; 
fltrl=vBasis(y ,  grdl,  ordl,  devl,  5,  grdlLength); 
fltr2=vBasis(y ,  grd2 ,  ord2 ,  dev2 ,  5,  grd2Length) ; 
wt(y,  5,  alpha); 
for (i=Q ; icnint ; i++) 

{ 

for(k=fltrl . index[Q] [i] ;  k<=fltrl . index[l] [i] ;  k++) 

{ 

for(j=fltr2.index[0] [i] ; j<=fltr2 . index [1] [i] ; j++) 

{ 

u=fltrl . v[k-fltrl . index[8] [i] ] [i] *fltr2 . v[ j -fltr2 . index[8] [i] ] [i] ; 
A [k] [j]=A[k] [j]+(b-a)*u*weight[i]*alpha[i]/2 ; 

} 

} 

} 

} 

return  A; 

} 

A.1.7  WEIGHTS 

#include  "spline. h" 

void  wt(double  y[] ,  int  length,  double  alpha[]) 

{ 


*  ©Author:  Dane  Burrows 
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*  ©Date  9-July-07 

*  ©Description: 

*  Returns  a  weight  for  calculating  the  inner  product. 


*  ©Usage: 

*  wt(y,  length,  alpha); 

*  Input : 

*  y  :  array  of  dimension:  length  which  can  be  used  for 

*  evaluating  the  weight. 

*  length  :  length  of  the  array:  y. 

*  alpha  :  empty  array  to  be  filled  with  the  result. 


forCint  i=0;i<5;i++) 
alpha[i]=l ; 

} 

A.1.8  PLOT  BASIS 

#include  "spline. h" 

Array2D<double>  plotBasis(Array2D<double>  v,  Array2D<int>  index,  double  x[] , 
int  order,  int  vLength,  int  ndim) 

{ 


*  ©Name:  Dane  Burrows 

*  ©Date  9-July-07 

*  ©Description: 

*  This  routine  evaluates  the  values  of  the  B-spline  basis. 

*  ©Usage: 

*  output  <name>=vBasis(x,  grid,  order,  derivative,  xLength,  gridLength) ; 

*  Input : 

*  v  :  an  array  of  dimension  order  +1  by  M,  where  M  is  the 

*  length 

*  of  the  array  x 

*  index  :  indices  of  the  basis  elements  with  non-zero  values 

*  at  a  point  x.  index  is  a  2  by  M  array, 

*  index[0][k]  --  lowest  index  of  non-zero  basis 

*  element  at  x[k] . 

*  index[l][k]  --  highest  index  of  non-zero  basis 

*  element  at  x[k] . 
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X 


:  array  of  values  for  x  on  which  the  basis  functions 
are  to  be  evaluated, 
order  :  order  of  the  spline  functions. 

vLength  :  length  of  the  array  v. 

ndim  :  total  number  of  basis  elements,  ndim=N+order-l . 

Output : 

u  :  an  array  that  can  be  graphed  to  show  the  basis 

functions . 


.  j.  j.  j.  .i.  j 


int  i,  j,  k; 

Array2D<double>  u(ndim,  vLength); 
for(i=8; i<vLength; i++) 

{ 

for ( j  ; j cndim ; j  ++) 

{ 

if(index[8] [i]<=j&&j<=index[l] [i]) 

{ 

u[j] [i]=v[j-index[8] [i] ] [i] ; 

} 

else 

u[j][i]=8; 

} 

} 

return  u; 

} 

A.1.9  SIGN 

#include  "spline. h" 

double  sign(double  a) 

{ 


*  ©Author:  Dane  Burrows 

*  ©Date  9-July-87 

*  ©Description: 

*  Checks  the  sign  of  a  number  and  returns  the  sign  as  either  1,  -1,  or  8. 

*  ©Usage: 

*  output  <name>=vBasis(x,  grid,  order,  derivative,  xLength,  gridLength) ; 

*  Input : 

*  a  :  A  double  value  to  have  its  sign  evaluated. 

*  Output : 

*  x  :  A  double  value  of  either  1.8  (positive)  -1.8 
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(negative)  or  ®.Q. 


double  x=8; 
if (a==Q) 
x=® ; 
if (a<8) 
x=-l; 
if (a>8) 
x=l; 

return  x; 


} 

A.1.10  SURF  VALUE 

#include  "spline. h" 

Array2D<double>  surf_value (Array lD<double>  x,  Array lD<double>  y, 

Array lD<double>  zgrid,  ArraylD<double>  nrgrid,  Array2D<double>  rgrid, 
int  rorder,  int  zdev,  int  rdev,  Array2D<double>  coefz) 

{ 

*  ©Author:  Dane  Burrows 

*  ©Date  9-July-®7 

*  ©Description: 

*  Evaluate  the  approximation  for  a  finite  number  of  slices  at  given  points 

*  zgrid[i]  using  interpolation  in  z. 

*  ©Usage: 

*  Array2D  <name>=surf_val(x,  y,  zgrid,  nrgrid,  rgrid,  rorder,  zdev,  rdev, 

*  coefz) ; 

*  Input : 

*  x  :  z  values  where  surface  is  requested. 

*  y  :  r  values  where  surface  is  requested. 

*  zgrid  :  grid  points  in  z. 

*  nrgrid  :  number  of  r  points  at  each  z  grid. 

*  rgrid  :  r  grid  points. 

*  order  :  order  of  spline. 

*  dev  :  order  of  derivative. 

*  coefz  :  coefficents  with  respect  to  the  standard  basis  at  each 

*  z  values  in  zgrid[i] . 

*  Output : 

*  v  :  value  of  the  z  interpolating  slice  function. 
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int  i,  j,  k,  lenx,  n; 
double  weightl,  weight2; 

Array2D<double>  v; 

n=zgrid. dim() ; 
lenx=x.dim()-l; 

for(i=®;  i<lenx;  i++) 

{ 

Array lD<double>  tmp; 
for(j=® ; j  <n- 1 ;  j++) 

{ 

//Interpolate  between  the  first  two  slices 

// 

Array lD<double>  tmprgridl((int)nrgrid[j]) , 
tmprgrid2((int)nrgrid[j+l]) , 
tmp_gridl((int)nrgrid[j] ) ,  tmp_grid2 ((int)nrgrid[j+l] ) , 
tmp_grid_z((int)nrgrid[j]) ,  tmp_grid_r((int)nrgrid[j] ) , 
tmp_coefl(coefz.dira2  0)  ,  tmp_coef2(coefz.dim2Q)  , 
tmp_coef_z(coefz.dim2Q)  ,  trap_coef_r (coefz . dim2 ())  ; 
for(k=®;k<nrgrid[j] ;k++) 

{ 

tmprgridl [k]=rgrid[j] [k] ; 
tmprgrid2 [k]=rgrid[j  +  l] [k]  ; 

} 

if((x[i]  >=  zgrid[j])&&(x[i]<zgrid[j+l])) 

{ 

weight l=(x[i] -zgrid[j] )/(zgrid[j+l] -zgrid[j] ) ; 
tmp_gridl=tmprgridl ; 
tmp_grid2=tmprgrid2 ; 
if (zdev==l) 

{ 

weight2=l/ (x[i]*(log(zgrid[j+l])-log(zgrid[j]))) ; 
for (k=® ; k<tmp_gridl . dimQ ; k++) 

tmp_grid_z[k]=(tmp_grid2 [k] -tmp_gridl [k] )*weight2 ; 

} 

for (k=®;k<tmp_gridl. dimQ ;k++) 

{ 

tmp_grid_r [k] =tmp_gridl [k]  +  (tmp_grid2 [k] -tmp_gridl [k] )  -'weight 1 

} 

for (k=® ; k<coefz . dim2 () ; k++) 

{ 

tmp_coefl[k]=coefz[j] [k] ; 
tmp_coef2 [k] =coefz [ j+1] [k] ; 

} 

if (zdev==l) 

{ 

weight2=l/ (x[i]*(log(zgrid[j+l])-log(zgrid[j]))) ; 


41 

DISTRIBUTION  A:  Approved  for  Public  Release ;  Distribution  Unlimited 


for(k=0;k<tmp_gridl.dim() ;k++) 

tmp_coef_z[k]=(tmp_coef2  [k]  -tmp_coef  1  [k] ) -'weight 2  ; 

} 

for(k=@;k<tmp_coefl .dimQ ;k++) 

{ 

tmp_coef_r [k]  =  tmp_coefl [k]+(tmp_coef2 [k] -tmp_coef 1 [k] )*weightl ; 

} 

tmp=vspline(y ,  tmp_grid_r,  rorder,  zdev,  tmp_coef_r) ; 

} 

} 

//  Calculate  the  surface  value  at  points  requested 

// 

if (i==Q) 

v=Array2D<double>  (tmp.dimO,  lenx) ; 
for(j=0;  jctmp.dimO  ;  j++) 

{ 

v[j] [i]=tmp[j] ; 

} 


return  v; 

} 

A.2  An  example  in  C 

A.2.1  VBASIS 

#include  "spline. h" 

int  main() 

{ 

int  order,  derivative; 
int  xlength,  i,  gridLength,  j; 
gridLength=ll ; 
xlength=101 ; 

Array lD<double>  grid(gridLength) ,  x(xlength) ; 

Array2D<double>  plots; 
ofstream  outputl("vBasis.out") ; 

for (i=0 ; i<gridLength ; i++) 
grid[i]=i ; 

for(i=0;  i<xlength;  i++) 
x[i]=i/10.0; 

order=3 ; 
derivative^©; 

output  vbout=vBasis(x,  grid,  order,  derivative,  xlength,  gridLength); 
plots=plotBasis(vbout . v,  vbout. index,  x,  order,  vbout.v.dim2() ,  vbout.ndim); 
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for ( j =® ; j  <vbout . ndim ; j  ++) 

{ 

for(i=Q;i<vbout.v.dim2() ;i++) 

{ 

outputl  «x[i]  «  "  "  «  plots[j]  [i]«endl; 

} 

outputl  «endl; 

} 

return  ®; 

} 

A.2.2  VBNEUMANN 

#include  "spline. h" 

int  main() 

{ 

int  order,  derivative; 
int  xlength,  i,  gridLength,  j; 
gridLength=ll ; 
xlength=l®l ; 

Array lD<double>  grid(gridLength) ,  x(xlength) ; 
for (i=® ; i<gridLength ; i++) 
grid[i]=i ; 

for(i=®;  i<xlength;  i++) 
x[i]=i/l®.8; 

ofstream  outputl("vbneumann. out") ; 

order=3 ; 
derivative^®; 

output  vbout=vbneumann(x,  grid,  order,  derivative,  xlength,  gridLength); 
Array2D<double>  plots  =  plotBasis(vbout . v,  vbout. index,  x,  order, 
vbout.v.dim2() ,  vbout. ndim); 

for (i=® ; i<vbout . ndim ; i++) 

{ 

for (j  =Q ;  j  cvbout . v . dim2 () ; j  ++) 

{ 

outputl  «  x [ j ]  «  "  "  «  plots[i][j]  «  endl; 

} 

outputl  «  endl; 

} 

return  ®; 

} 

A.2.3  V  SPLINE 

#include  "spline. h" 
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int  main() 

{ 

int  order,  derivative; 

int  xlength,  i,  xdataLength,  gridLength,  zLength=ll,  yLength=101, 
aLength=2 1 ,  j ; 
xdataLength=ll ; 
gridLength=ll ; 
xlength=l®l ; 

ofstream  outputl("ydata.out") ; 
ofstream  output2 ("vspline . out") ; 

Array lD<double>  xdata(xdataLength) ,  grid(gridLength) ,  ydata(xdataLength) , 
x(xlength),  vspln; 

for (i=® ; i<xdataLength ; i++) 
xdata[i]=i ; 

for (i=® ; i<gridLength ; i++) 
grid[i]=i ; 

for(i=®;  i<xdataLength;  i++) 

{ 

ydata[i]=exp(-xdata[i]) ; 

outputl  «  xdata[i]  «  "  "  «  ydata[i]  «  endl ; 

} 

for(i=®;  i<xlength;  i++) 
x[i]=i/l®.®; 
order=3 ; 
derivative^®; 

Array lD<double>  approx; 

approx=lsqapp(xdata,  ydata,  ydata,  grid,  order,  8.1,  ®.®1); 
vspln=vspline(x,  grid,  order,  derivative,  approx); 

for(int  i=Q;  i  <  vspln. dim();  i  ++) 

output2  «  i/((double)xlength-l)*10.®  «  "  "  «  vspln[i]«  endl; 

return  ®; 

} 

A.2.4  LEAST  SQUARES  APPROXIMATION 

#include  "spline. h" 

int  main() 

{ 

int  order; 

int  xlength,  i,  j,  xdataLength,  gridLength; 
xdataLength=ll ; 
gridLength=ll ; 

Array lD<double>  xdata(xdataLength) ,  grid(gridLength) ,  ydata(xdataLength) ; 
for (i=® ; i<xdataLength ; i++) 
xdata[i]=i ; 

for (i=® ; i<gridLength ; i++) 
grid[i]=i ; 

for(i=8;  i<xdataLength;  i++) 
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ydata[i]=exp(-xdata[i]) ; 
ofstream  outputl("lsqapp.out") ; 
order=3 ; 

Array lD<double>  approx; 

approx=lsqapp(xdata,  ydata,  ydata,  grid,  order,  8.1,  Q.®1); 

for(i=0; i<approx.dim() ; i++) 

{ 

outputl  «  approx[i]  «  endl; 

} 

return  0; 

} 

A.2.5  SLICE  COEFFICIENTS 

#include  "spline. h" 

int  main() 

{ 

int  order,  derivative; 

int  i,  zLength=ll,  aLength=21,  j; 

order=3 ; 
derivative^©; 

Array lD<double>  z(zLength),  nr(zLength) ; 
for(i=0;  i<zLength;  i++) 
nr [i]=aLength; 

Array2D<double>  r(zLength,  aLength) ; 

Array2D<double>  nval(zLength,  aLength); 
ofstream  outputl("slice.out") ; 

for(i=0;  i<zLength;  i++) 
z[i]=i; 

for(i=0;  i<zLength;  i++) 

{ 

for (int  j=0;  j<nr[0] ;  j++) 

{ 

r[i] [j]=-10+j ; 

nval [i] [j]=-0. 001*exp(-r [i] [j]*r [i] [j])*exp(-z [i] ) +1 .35; 

} 

} 

Array2D<double>  slice; 

slice=slice_coef(z ,  nr,  r,  nval,  z,  nr,  r,  order,  .1,  .01); 
for(i=0;i<slice.dim2() ; i++) 

{ 

for(j=0; j<slice.diml() ; j++) 

{ 

outputl  «  j  «  "  "  «  i  «  "  "  «  slice[j][i]  «  endl; 

} 
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output 1  «  endl; 


} 

return  ®; 

} 

A.2.6  SURF  VALUE 

#include  "spline. h" 

int  main() 

{ 

int  order,  derivative; 

int  xlength,  i,  gridLength,  zLength=ll,  yLength=l®l,  aLength=21,  j; 
xlength=l®l ; 

Array lD<double>  x(xlength),  y(yLength) ; 
for(i=®;  i<xlength;  i++) 
x[i]=i/l®.®; 

for(i=®;  i<yLength;  i++) 
y[i]=-l®+i/5.8; 
ofstream  outputl("surf . out") ; 

order=3 ; 
derivative^®; 
double  temp[l]; 

Array lD<double>  z(zLength),  nr(zLength) ; 
for(i=®;  i<zLength;  i++) 
nr [i]=aLength; 

Array2D<double>  r(zLength,  aLength); 

Array2D<double>  nval(zLength,  aLength); 

for(i=®;  i<zLength;  i++) 
z [i]=i ; 

for(i=®;  i<zLength;  i++) 

{ 

for (int  j=Q;  j<nr[Q] ;  j++) 

{ 

r[i] [j]=-l®+j ; 

nval [i] [j]=-® . 8®l*exp(-r [i] [j]*r [i] [j])*exp(-z [i])+l . 35 ; 

} 

} 

Array2D<double>  slice,  surf; 

slice=slice_coef(z ,  nr,  r,  nval,  z,  nr,  r,  order,  .1,  .01); 
surf=surf_value(x,  y,  z,  nr,  r,  order,  ®,  ®,  slice); 
for(i=Q; i<surf .dim2() ;i++) 

{ 

for(j=®; j<surf .diml() ; j++) 

{ 

outputl  «  j  «  "  "  «  i  «  "  "  «  surf[j][i]  «  endl; 

} 
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output 1  «  endl; 


} 

return  ®; 

} 

A.3  MATLAB  Code 

A.3.1  VBASIS 

function  [v , ndim , index] =vbasis (x , grid , order , derivative) 

% 

%  Description:  This  routine  evaluates  the  values  of  the  B-spline  basis 
%  (or  the  derivatives)  functions  at  given  points.  The  grid  points  and 

%  the  order  of  the  splines  are  specified  by  user,  however,  additional 

%  grid  points  outside  of  the  interval  [xmin,  xmax]  are  choosen  by 
%  the  program  to  provide  a  complete  basis. 

%  Usage: 

%  [v, ndim, index] =vbasis(x, grid, order  derivative) ; 

%  Input : 

%  x  : 

% 

%  grid  : 

% 

% 

% 

% 

%  order  : 

%  derivative  : 

% 

%  Output : 

%  v  : 

% 

%  ndim  : 

%  index  : 

% 

% 

% 

% 

% 

% 

N=length(grid) ; 

M=length(x) ; 
ndim=N+order-l ; 

% 

%  Construct  the  augmented  grid. 

% 

dgrid=grid(2) -grid(l) ; 

lgrid=grid( 1) -order*dgrid : dgrid : grid(l) -dgrid ; 
dgrid=grid(N)-grid(N-l) ; 

rgrid=grid(N)+dgrid : dgrid : grid(N) +order "dgrid ; 


row  vector  of  values  for  x  on  which  the  basis  functions 
are  to  be  evaluated. 

the  grid  points  in  ascending  order,  all  grid  points 
must  be  distinct.  The  interval  on  which  the  spline 
basis  functions  are  defined  are  given  by: 

[grid(l) ,grid(N)] 

where  N  is  the  length  of  the  vector  grid, 
order  of  the  spline  functions, 
order  of  derivative  needed. 

derivative=®,  value  of  the  function  is  requested. 

an  array  of  dimension  order+1  by  M,  where  M  is  the 
length  of  vector  x. 

total  number  of  basis  elements,  ndim  =  N+order-1. 
the  indices  of  basis  elements  with  non-zero  values 
at  a  point  x.  index  is  a  2  by  M  matrix, 

index (l,k)  --  lowest  index  of  non-zero  basis  element 

at  x(k) . 

index(2,k)  --  highest  index  of  non-zero  basis  element 
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agrid=[lgrid,  grid,  rgrid] ; 


Main  loop  over  points  x. 

v=[] ; 
index=[] ; 
for  k=l :M 

Determine  the  interval  [grid(j),  grid(j+l)]  in  which  x(k)  belongs, 
for  j=l:N-l 

if  sign(x(k)-grid(j))*sign(grid(j+l)-x(k))  >=  ® 
break; 
end; 
end; 

if  x(k)  <  grid(l) 

j=i ; 

end; 

if  x(k)  >  grid(N) 
j=N- 1 ; 
end; 

Evaluate  the  values  of  the  basis  functions  (or  derivatives)  at  x(k) . 

1.  Evaluate  the  values  of  the  truncated  polynomials. 

factor=l ; 
if  derivative  >  ® 

for  i=0: derivative- 1 

factor=factor*(order-i) ; 

end; 

end; 

trunc  =  zeros(l , 2*order+2) ; 
localg  =  agrid(j : j+2*order+l) ; 
for  i=l : 2*order+2 

trunc (i)=max( [(x(k)-agrid(j+i-l)) , ®] ) ; 
if  order  >  derivative 

trunc(i)=factor*(trunc(i) ~ (order-derivative)) ; 
elseif  order  ==  derivative 

trunc (i) =f actor *sign (trunc (i)) ; 
else 

disp(’The  spline  function  is  not  differentiable’); 
return; 
end; 

end; 

2.  Compute  the  divided  differences. 

for  l=l:order+l 

for  11=1 : 2*order+2-l 

trunc (11)= (trunc (11+1) -trunc (11)) /(localg (11+1) -localg (11)) ; 

end; 
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end; 

% 

%  3.  Store  the  value  in  the  vector  v. 

% 

itrunc=[j ; j+order] ; 
index  =  [index  itrunc] ; 
trunc  =  trunc(l :order+l) ; 
v  =  [v  trunc ’ ] ; 
end; 

A.3.2  VBNEUMANN 

function  [v , ndim , index] =vbneumann(x , grid , order , derivative) 

% 

%  Description:  evaluates  the  value  of  basis  elements  of  spline  functions 
%  of  the  specified  order  on  the  given  grid  which  satisfies  the 
%  Neumann  boundary  conditions. 

%  Usage: 

%  [v , ndim , index] =vbneumann(x , grid , order , derivative) ; 

%  Input : 

%  x  : 

% 

%  grid  : 

% 

% 

% 

% 

%  order  : 

%  derivative  : 

% 

%  Output : 

%  v  : 

% 

%  ndim  : 

%  index  : 

% 

% 

% 

% 

% 

% 

N  =  length(grid) ; 
interval  =  [grid(l) , grid(N) ] ; 

[tau, ndim, index]  =  vbasis(interval , grid, order , 1) ; 

[u, ndim, index]  =  vbasis (x, grid, order, derivative) ; 

M  =  length(x) ; 
v=u; 

% 

%  Impose  the  boundary  conditions. 

% 

for  k=l:M 


row  vector  of  values  for  x  on  which  the  basis  functions 
are  to  be  evaluated. 

the  grid  points  in  ascending  order,  all  grid  points 
must  be  distinct.  The  interval  on  which  the  spline 
basis  functions  are  defined  are  given  by: 

[grid(l) ,grid(N)] 

where  N  is  the  length  of  the  vector  grid, 
order  of  the  spline  functions, 
order  of  derivative  needed. 

derivative=®,  value  of  the  function  is  requested. 

an  array  of  dimension  order+1  by  M,  where  M  is  the 
length  of  vector  x. 

total  number  of  basis  elements,  ndim  =  N+order-3. 
the  indices  of  basis  elements  with  non-zero  values 
at  a  point  x.  index  is  a  2  by  M  matrix, 

index Cl,k)  --  lowest  index  of  non-zero  basis  element 
at  x(k) . 

index(2,k)  --  highest  index  of  non-zero  basis  element 
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if  index(l,k)  ~=  1 

index(l,k)  =  index(l ,k)-l ; 
else 

for  i=l: order 

v(i,k)  =  u(i+l,k)-tau(i+l, l)*u(l ,k)/tau(l , 1) ; 

end; 

end; 

if  index(2,k)  "=  ndira 

index(2,k)  =  index(2 ,k)-l ; 
else 

for  i=l: order 

v(i,k)  =  v(i ,k)-tau(i , 2)*u(order+l ,k)/tau(order+l , 2) ; 

end; 

index(2,k)  =  ndim-2; 
end; 

end; 

ndim  =  ndira- 2; 
return; 

A.3.3  LEAST  SQUARES  APPROXIMATION 

function  [coef]=lsqapp(xdata,ydata,wdata, . . . 

xgrid , order , filter , alpha® , alphal) 

% 

%  Description:  Compute  the  least  square  approximation  of  the  data  set 
%  using  a  given  set  polynomial  spline  functions.  The  optimization 
%  functional  is  given  by 

%  J(coef)  =  \sum"N_{ j=l}  wdata_j | L(t_j)-S(t_j) | "2 
%  \alpha_®\int"{t_max}_{t_min}  | L(t) -S(t) | "2dt 

%  \alpha_l\int"{t_max}_{t_min}  |L’ (t)-S’ (t) | "2dt, 

%  where : 

%  N  :  number  of  data  points. 

%  L  :  linear  spline  interpolation  of  the  data. 

%  S  :  polynomial  spline  function. 

%  \alpha_®  :  weight  on  the  L_2  norm. 

%  \alpha_l  :  weight  on  the  H_1  norm. 

%  Usage: 

%  [coef ] =lsqapp (xdata , ydata , wdata , xgrid , order .filter , alpha® , alphal) 

%  where : 

%  xdata  :  data  values  for  the  independent  variable. 

%  ydata  :  data  values  for  the  dependent  variable. 

%  wdata  :  weights  on  the  data  points. 

%  xgrid  :  grid  for  the  spline  function. 

%  order  :  order  of  the  polynomial  spline. 

%  filter  :  routine  for  evaluation  of  the  spline  functions. 

%  alpha®  :  weight  on  the  L_2  norm. 

%  alphal  :  weight  on  the  H_1  norm. 

%  coef  :  coefficients  for  the  optimal  spline  function. 

% 

[u.ndata, index]=vbasis(xdata, xdata, 1,Q) ; 

% 
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%  Evaluate  the  pointwise  term. 

% 

PI  =  zeros (ndata,ndata) ; 
for  k=l:ndata 

for  j=index(l,k) :index(2,k) 

Pl(k,j)  =  u(j-index(l ,k)+l ,k) ; 

end; 

end; 

intp  =  inv(Pl)*ydata’ ; 

[u, ndim , index] =feval (filter , xdata , xgrid , order , ®) ; 

Q  =  zeros (ndata, ndim) ; 
for  k=l: ndata 

for  j=index(l,k) :index(2,k) 

Q(k, j)  =  u(j -index (1 , k)+l,k) ; 

end; 

end; 

W  =  diag(wdata, ®) ; 
rl  =  Q’*W*ydata’; 

A1  =  Q ’ *W*Q ; 

% 

%  Evaluate  the  L_2  term. 

% 

xmin  =  min( [xdata  xgrid]); 
xmax  =  max([xdata  xgrid]); 

[P2]  =  innprd(xgrid, order, ®, filter, xdata, 1,8, 'vbasis’ ,.. . 
’unif ’ , xmin, xmax) ; 

r2  =  P2*intp; 

[A2]  =  innprd(xgrid, order,®, filter, xgrid, order ,.. . 

®, filter, ’unif’ , xmin, xmax) ; 

% 

%  Evaluate  the  H_1  term. 

% 

[P3]  =  innprd(xgrid, order , 1 , filter , xdata, 1 , 1 , 'vbasis ’,.. . 
’unif’ , xmin, xmax) ; 

r3  =  P3*intp; 

[A3]  =  innprd(xgrid, order, 1, filter, xgrid, order ,.. . 

1, filter, ’unif’ , xmin, xmax) ; 

% 

%  Solve  for  the  optimal  coefficients. 

% 

r  =  rl+alpha®*r2+alphal*r3 ; 

A  =  Al+alpha®*A2+alphal*A3 ; 
coef  =  (inv(A)*r)’; 
return; 

A.3.4  V  SPLINE 

function  [v] =vspline (x , xgrid , order , dev , filter , coef) 

% 

%  Description:  Evaluate  a  given  polynomial  spline  function. 

%  Usage: 
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%  [v]=vspline(x, xgrid, order, dev, filter, coef) ; 


%  where: 

%  x 

%  xgrid 

%  order 

%  dev 

%  filter 

%  coef 

%  v 


values  of  the  independent  variable, 
grid  on  which  the  splines  are  defined, 
order  of  spline, 
order  of  derivative. 

filter  used  to  impose  boundary  conditions, 
coefficients  with  respect  to  the  standard  basis, 
value  of  the  spline. 


% 


[u , ndim , index] =f eval (filter , x , xgrid , order , dev) ; 
if  ndim  ~=  length(coef) 

dispC’The  dimension  of  the  coefficient  vector  is  wrong.’); 
return; 
end; 

v  =  zeros(size(x)) ; 


%  Calculate  the  spline  values. 

% 

for  k=l : length(x) 

for  j=index(l,k) :index(2,k) 

v(k)  =  v(k)+coef(j)*u(j-index(l,k)+l,k) ; 

end; 

end; 
return ; 


A.3.5  2D-SLICE  COEFFICIENT 

function  [coefp]=slice_coef (p ,nt , t , ival ,pgrid,ntgrid, tgrid, order , filter , . . 

alpha®, alphal) 

% 

%  Description:  Evaluate  coefficient  for  each  slice  at  p_{ i }  to  approximate 
%  the  data  set  t(i,l:nt(i))  and  ival(i , 1 :nt(i)) .  This  approximation 
%  is  done  using  one  dimensional  approximation. 

%  Usage: 

%  [coef]  =  slice_coef(pgrid, ntgrid, tgrid, order, filter, p,nt,t, ival, . . 

%  alpha®, alphal) ; 

% 

%  Inputs: 

%  pgrid 

%  ntgrid 

%  tgrid 

%  order 

%  filter 

%  p 

%  nt 

%  t 

%  ival 

%  alpha® 

%  alphal 

% 


grid  points  in  z 

number  of  grid  points  in  r 

grid  points  in  r 

order  of  the  spline  requested. 

filter  requested. 

z  values  for  data. 

number  r  data  points  at  each  slice. 

r  measurements. 

refractive  index  measurements 

weight  for  data  approximation. 

weight  for  derivative  approximation. 
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:  lenp  by  (nt(i)+order- 1)  matrix  containing  the 
approximation  coefficients  at  each  slice. 


%  Outputs: 

%  coefp 

% 

% 

lenp  =  max(size(p)) ; 
coefp  =  [] ; 
for  i=l:lenp 

[tmp]=lsqapp(t(i , l:nt(i)) ,ival(i, l:nt(i)) ,ones(l ,nt(i)) , . . . 

tgrid(i, l:ntgrid(i)) , order , filter , alphas , alphal) ; 
coefp  =  [coefp ;tmp]; 

end; 
return ; 

A.3.6  INNER  PRODUCT 


function  [A]=innprd(grdl , ordl , devl , fltrl , grd2 , ord2 , dev2 , fltr2 , wt , xmin , xmax) 

% 

%  Description:  Compute  the  matrix  of  the  inner  product  of  two  families  of 
%  polynomial  spline  basis  functions.  If  {B_k} ,  k=l,...,N  and  {C_j}, 

%  j=l . M,  then  the  matrix  is  given  by 

%  A_{k,j}  =  <B_k,  C_j>. 

%  Usage: 

%  [A]=innprd(grdl,ordl,devl, ’fltrl’ , grd2 , ord2 , dev2 , ’ fltr2 ’ , wt , xmin , xmax) ; 
%  Input : 

%  grdl 

% 

%  ordl 

%  devl 

% 

%  fltrl 

% 

%  grd2 

% 

%  ord2 

%  dev2 

% 

%  fltr2 

% 

%  wt 

%  xmin 

%  xmax 

%  Output : 

%  A 

% 

weight  =  [8.2369268851,  8.4786286785,  8.5688888889]; 
weight  =  [weight  weight(2)  weight(l)] ; 
x  =  [-8.9861798459,  -8.5384693181,  8]; 
x  —  [x  -x(2)  -x(l)] ; 

% 

%  Construct  the  combined  grid. 


:  grid  points  of  the  first  group  of  polynomial  spline 
functions . 

:  the  order  of  the  first  group  of  spline  functions. 

:  order  of  the  derivatives  of  the  first  group  of 
spline  functions. 

:  filter  that  modifies  the  basis  elements  in  order 
to  satisfy  the  boundary  condition. 

:  grid  points  of  the  second  group  of  polynomial  spline 
functions . 

:  the  order  of  the  second  group  of  spline  functions. 

:  order  of  the  derivatives  of  the  second  group  of 
spline  functions. 

:  filter  that  modifies  the  basis  elements  in  order 
to  satisfy  the  boundary  condition. 

:  weight  function  in  the  inner  product. 

:  lower  bound  of  the  interval  of  integration. 

:  upper  bound  of  the  interval  of  integration. 

:  the  matrix  of  inner  product. 


53 

DISTRIBUTION  A:  Approved  for  Public  Release ;  Distribution  Unlimited 


s?  s?  s? 


% 


nint  =  5 ; 
n  =  length(grdl) ; 
dgrd=grdl (2) -grdl (1) ; 

lgrd=grdl ( 1) -ordl*dgrd : dgrd : grdl (1) -dgrd ; 
dgrd=grdl(n)-grdl(n-l) ; 

rgrd=grdl (n)+dgrd : dgrd : grdl (n) +ordl*dgrd ; 
agrdl=[lgrd,  grdl,  rgrd] ; 
n  =  length(grd2) ; 
dgrd=grd2 (2) -grd2 (1) ; 

lgrd=grd2 (1) -ord2*dgrd : dgrd : grd2 (1) -dgrd ; 
dgrd=grd2 (n) -grd2 (n- 1) ; 

rgrd=grd2 (n)+dgrd : dgrd : grd2 (n) +ord2*dgrd ; 
agrd2=[lgrd,  grd2 ,  rgrd]; 
cgrd  =  [agrdl  agrd2] ; 
cgrd  =  sort(cgrd); 

% 

%  Choose  the  maximum  and  the  minimum  grid  point. 

% 

igrd  =  xmin; 
k  =  1; 

for  cx  =  cgrd 

if  cx  >=  xmax 
break; 
end; 

if  igrd(k)  <  cx 

igrd  =  [igrd  cx] ; 
k  =  k+1; 
end; 

end; 

igrd  =  [igrd  xmax] ; 

[vl,Nl,indexl]  =  feval(fltrl , grdl (1) , grdl , ordl , devl) ; 
[v2,N2,index2]  =  feval(fltr2 , grd2 (1) , grd2 , ord2 , dev2) ; 

Calculate  the  inner  product  matrix. 

A  =  zeros(Nl ,N2) ; 
for  k=  l:length(igrd)-l 
a  =  igrd(k) ; 
b  =  igrd(k+l) ; 
y  =  (b+a)/2+(b-a)*x/2 ; 

[vl ,ndim, indexl]  =  feval(fltrl , y , grdl , ordl , devl) ; 
[v2 ,ndim, index2]  =  feval(fltr2 , y , grd2 , ord2 , dev2) ; 
alpha  =  feval(wt.y); 
for  i  =  l:nint 

for  k  =  indexl (1, i) : indexl (2, i) 
for  j  =  index2(l,i) :index2(2,i) 

u  =  vl(k-indexl(l ,i)+l ,i)*v2(j-index2(l , i)+l , i) ; 
A(k,j)  =  A(k, j)+(b-a)*u*weight(i)*alpha(i)/2 ; 
end; 
end; 
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end; 

end; 

return; 

A.3.7  UNIFORM  WEIGHTS 

function  [v]=unif(x) 

% 

%  Description: 

%  Uniform  weight  function. 

% 

v  =  ones(size(x)) ; 
return; 

A.3.8  PLOT  BASIS 

function  plot_basis(x , knots , order , deriv , filter) 

%  This  function  plots  the  basis  functions  for  a  particular  polynomial  splines 
%  with  the  given  order  and  derivative  for  the  knots  specified. 

[v,ndim, index]  =  feval(filter ,x,  knots,  order,  deriv); 
nx  =  length(x) ; 
v2  =  zeros (ndim, nx) ; 
for  j=l:nx 

for  i=index(l , j) : index (2 , j) 

v2(i,j)  =  v (i- index (1 , j)+l , j) ; 

end; 

end; 

plot(x , v2 (1 , :)) ; 
hold  on; 
pause ; 

for  i=2:ndim 

plot(x,v2(i, :)) ; 
pause; 

end; 

return; 


A.3.9  SURF  VALUE 


function  [v] =surf_value (x , y , pgr id , ntgrid , tgrid , torder , pdev , tdev , filter , coefp) 

% 


%  Description:  Evaluate  the  approximation  for  a  finite  number  of  slices  at 
%  given  z  points  pgrid(i)  using  logarithmic  interpolation  in  p. 

%  Usage: 

%  [v]  =  surf_val (x,y,pgrid, ntgrid, tgrid, order, dev, filter, coefp) ; 


%  where: 

%  x 

%  y 

%  pgrid 

%  ntgrid 

%  tgrid 

%  order 

%  dev 


z  values  where  surface  is  requested, 
r  values  where  surface  is  requested, 
grid  points  in  z. 

number  of  r  grid  points  at  each  z  grid, 
r  grid  points, 
order  of  spline, 
order  of  derivative. 
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%  filter  :  filter  used  to  impose  boundary  conditions. 

%  coefp  :  coefficients  with  respect  to  the  standard  basis  at  each 

%  pressure  values  in  pgrid(i) . 

%  v  :  value  of  the  p-logarithmic  interpolating  slice  function. 

% 

lenx  =  length(x) ; 
v  =  []; 

n  =  length(pgrid) ; 
for  i=l:lenx 
for  j=l:n-l 

% 

%  Interpolate  between  first  two  slices. 

% 

if  ((x(i)  >=pgrid(j))  &  (x(i)  <  pgridCj+1))) 

%weightl=(log(x(i))-log(pgrid(j)))/(log(pgrid(j+l))- 
log(pgridCj))) ; 

weight l=(x(i)-pgrid(j))/ (pgrid(j+l)-pgrid(j)) ; 
tmp_gridl  =  tgridCj , 1 :ntgrid(j)) ; 
tmp_grid2  =  tgrid(j+l , 1 :ntgrid(j+l)) ; 
if  (pdev  ==  1) 

weight2  =  l/(x(i)*(logCpgridCj+l))-logCpgrid(j)))) ; 
tmp_grid_p  =  (tmp_grid2-tmp_gridl)*weight2 ; 
end; 

tmp_grid_t  =  tmp_gridl  +  (tmp_grid2-tmp_gridl)  *weightl ; 
tmp_coefl  =  coefpCj,:); 
tmp_coef2  =  coefp(j+l , : ) ; 
if  (pdev  ==  1) 

weight2  =  l/(x(i)*(log(pgrid(j+l))-log(pgrid(j)))) ; 
tmp_coef_p  =  (tmp_coef2-tmp_coefl)*weight2 ; 
end; 

tmp_coef_t  =  tmp_coefl  +  (tmp_coef2-tmp_coefl)*weightl ; 
end; 
end; 

% 

%  Calculate  the  surface  value  at  points  requested. 

% 

tmp  =  vspline (y,tmp_grid_t,torder,tdev, filter, tmp_coef_t) ; 
v  =  [v , tmp ’ ] ; 

end; 

return; 

A.3.10  DEFINING  THE  REFRACTIVE  INDEX 

function  [n]  =  rindex_profile(r,z,mua,theta,A,nbar) 

%This  function  computes  the  refractive  index  profile  from  this  class  of 
functions: 

%  n(r,z)  =  A*exp(-r"2/theta''2)*exp(-mua*z) 

% 

%  Input: 

%  x  :  x  grid  of  size  nx  x  1 

%  y  :  y  grid  of  size  ny  x  1 
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z 


% 

% 

% 

% 

% 

% 


mu  a 

theta 

A 

nbar 


%  Output: 
%  n 


% 


z  grid  of  size  nz  x  1 
absorption  coefficient  (scalar) 
spread  parameter  (scalar) 
percent  deviation  from  nbar  (scalar) 
background  refractive  index  (scalar) 


refractive  index  profile  at  (r,z)  of  size  nr  x  nz 


nr  =  length(r) ; 
nz  =  length(z) ; 
n  =  zeros (nr ,nz) ; 
for  i=l:nr 

for  j=l:nz 

n(i , j)  =  -A*exp(-(r(i) "2)/(theta"2))*exp(-mua*z(j))+nbar ; 

end; 


end; 

return; 
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